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Preface 


The  efficient  calculation  of  flows  with  viscous/inviscid 
interaction  has  been  the  topic  of  much  research  over  the  past 
several  years.  Three  algorithms  were  developed  to  solve  the 
viscous/ inviscid  problem.  The  first  method  uses  finite 
difference  equations  with  successive  line  overfrelaxation 
(SLOP)  sweeps  for  solving  the  approximate  Navier-Stokes 
equations  in  the  viscous  region  and  the  stream  function 
equation  in  the  inviscid  region.  An  implicit  coupling  scheme 
is  developed  to  match  the  two  solutions.  The  second  method 
uses  finite  difference  approximations  for  solving  the  stream 
function  equation  in  the  inviscid  region  and  a  fourth  order 
Runge-Kutta  method  for  solving  the  integral  boundary  layer 
equations  in  the  viscous  region.  In  the  third  method,  the 
inviscid  flow  solution  is  obtained  by  a  panel  method,  while 
the  viscous  flow  solution  is  obtained  using  the  finite 
difference  form  for  the  boundary  layer  equations 
operating  in  an  inverse  scheme. 

I  would  like  to  thank  Dr.  A.  Halim  for  his  extensive 
technical  assistance  with  this  work  and  for  supervising  this 
thesis,  and  Dr,  Shang  of  the  Flight  Dynamics  Laboratory  for 
his  encouragement  and  insightful  comments.  Special  thanks 
goes  to  Lisa,  my  fiancee,  for  her  understanding  and  support 
over  the  last  year. 


The  ASD  Cyber  was  used  extensively  as  the  computational 
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(AFWAL/FIG)  for  sponsoring  the  work  performed  on  the  Cyber 
and  for  allowing  me  the  time  away  from  my  normal  duties  to 
conduct  this  research. 
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Abstract 

The  study  of  flows  with  viscous/inviscid  interaction  has 
attracted  many  researchers  over  the  last  decade.  These  flows 
occur  whenever  the  adverse  pressure  gradient  is  large  enough 
to  cause  flow  separation.  The  current  emphasis  is  to  find 
efficient  ways  of  solving  these  types  of  flows  without 
solving  the  full  Navier-Stokes  equations. 

Three  methods  for  solving  the  viscous/inviscid  problem 
were  studied.  The  first  method  uses  finite  difference 
equations  to  model  both  the  viscous  and  inviscid  regions.  A 
coupling  scheme  is  developed  to  match  the  two  solutions.  The 
second  method  solves  the  integral  boundary  layer  equations  in 
the  viscous  region  and  finite  difference  equations  in  the 
inviscid  region.  The  third  method  solves  the  Hilbert 
integral  to  generate  a  correction  to  the  inviscid  velocity 
using  the  boundary  layer  equations  as  the  viscous  model.  The 
model  problem  used  in  this  work  is  Howarth  flow  over  a  flat 
plate . 


The  three  methods  were  evaluated  in  terms  of  solution 
accuracy,  memory  requirements,  and  computation  times.  „ 
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NUMERICAL  STUDY  OF  THREE 
VISCOUS  /  INVISCID  INTERACTION  METHODS 

I  Introduction 

The  aerodynamic  performance  of  a  flight  vehicle  is 
greatly  affected  by  the  accuracy  of  the  computational  or 
experimental  methods  used  In  the  design  process.  Any 
successful  design  must  carefully  account  for  the  drag.  The 
VISCOUS  calculations  have  a  substantial  Impact  on  the  drag 
estimation.  Naturally,  the  full  Navler-Stokes  (NS)  equations 
correctly  estimate  the  drag  on  a  flight  vehicle.  However, 
their  use  can  be  costly  and  In  many  cases  unnecessary. 

Recent  research  has  focused  on  the  use  of  alternative 
forms  of  the  NS  equations.  These  approximate  sets  of 
equations  are  simpler  and  require  fewer  computations  than  the 
NS  equations,  but  are  valid  only  as  long  as  their  simplifying 
assumptions  are  not  violated.  For  example,  high  Reynolds 
number  flows  over  a  body  usually  result  in  the  formation  of  a 
thin  shear  layer  close  to  its  surface.  For  this  type  of  flow 
the  pressure  gradient  normal  to  the  body  and  viscous  terms 
with  derivatives  in  the  streamwise  direction  can  be  neglected 
by  an  order  of  magnitude  analysis  on  the  NS  equations.  The 
resulting  Boundary  Layer  (BL)  equations  are  widely  used  for 
many  high  Reynolds  number  flows.  One  of  the  attractive 
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features  of  the  BL  equations  is  that  they  are  parabolic. 

This  implies  that  the  solution  can  be  marched  in  the 
streair.wise  direction  without  iteration. 

Another  simplified  form  of  the  NS  equations  is  the 
Approximate  Navier-Stokes  (ANS)  equations.  The  ANS  equations 
assume  that  only  the  viscous  terms  with  derivatives  in  the 
streamwise  direction  are  small;  all  other  terms  are  retained. 
The  ANS  equations  fall  between  the  Navier-Stokes  and  the 
Boundary  Layer  equations  in  terms  of  accuracy.  They  are 
useful  because  they  are  partially  parabolized  in  the  case  of 
subsonic  flow  and  are  a  mixed  set  of  parabolic/hyperbolic 
equations  for  supersonic  flow  [1].  The  parabolized  equations 
allow  for  forward  marching  in  the  streamwise  direction.  For 
subsonic  flow,  forward  marching  is  still  possible,  however, 
several  iterations  may  be  necessary  to  achieve  convergence 
since  the  equations  still  contain  elliptic  inertia  type 
terms . 

Approximate  forms  of  the  NS  equations  are  particularly 
useful  when  used  in  recently  developed  zonal  techniques. 

These  techniques  divide  the  flow  region  into  distinct  zones, 
each  having  a  particular  set  of  assumptions  about  the  flow. 
For  example,  when  the  Reynolds  number  is  large  the  flow 
region  can  be  broken  up  into  a  viscous  region  (where  the  ANS 
or  BL  equations  can  be  used)  and  an  inviscid  region  where  a 
greatly  simplified  inviscid  model  is  used.  If  the  viscous 
region  is  small  in  comparison  with  the  inviscid  region  the 
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computational  savings  can  be  substantial.  This  is  true 
because  the  mathematical  model  for  the  viscous  flow  will  be 
solved  over  a  relatively  small  region  while  the  larger 
inviscid  region  uses  a  simple  inviscid  code.  A  coupling 
scheme  is  employed  in  the  zonal  technique  to  insure 
compatibility  between  the  two  regions.  The  iteration  of  the 
boundary  condition  at  the  interface  of  the  two  regions  Is  the 
mechanism  through  which  the  viscous  and  inviscid  flow  regions 
interact.  There  are  many  approaches  available  for  solving 
such  problems.  The  interacting  boundary  layer  theory  (IBLT) 
was  used  by  Carter  12],  Vatsa  et.  al.I3],  Edwards  and  Carter 
14]  and  Houwink  and  Veldman  [5]  to  solve  the  viscous/inviscid 
problem.  In  the  IBLT,  the  viscous  region  is  represented  by 
the  boundary  layer  equations;  the  inviscid  flow  can  be 
represented  in  a  number  of  different  ways  depending  upon  the 
flow  configuration  and  the  Mach  number.  Rubin  et.  al.l6]  and 
Swanson  et.  al.(7]  used  a  composite  velocity  representation 
of  the  inviscid  and  viscous  flow  regions.  Halim  and  Hafez 
[8-9]  solved  the  viscous  and  the  inviscid  regions  using  a 
semi-implicit  coupling  technique.  More  recent  efforts  were 
also  successful  using  a  fully  implicit  coupling  method  to 
obtain  efficient  solutions  [10-11 J  for  viscous/inviscid 
problems. 
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Any  numerical  method  developed  must  be  able  to  generate 
regular  solutions  in  the  event  of  an  adverse  pressure 
gradient.  If  separation  occurs  there  is  a  restriction  on 
what  can  be  used  as  the  boundary  condition  to  insure  regular 
behavior.  It  is  known  that  the  solution  of  the  boundary 
layer  equations  with  a  prescribed  pressure  gradient  results 
in  singularity  at  the  point  of  separation  [12J.  The 
singularity  at  the  separation  point  is  independent  of  the 
form  cf  the  equations  (i.e.  integral  or  differential).  The 
work  of  Carter  [13]  showed  that  regular  solutions  can  be 
obtained  with  an  inverse  approach  in  which  either  the 
displacement  thickness  or  the  skin  friction  is  specified.  In 
addition,  the  criterion  of  Meksyn  114]  states  that  if  (the 
velocity  at  the  outer  edge  of  the  viscous  layer)  did  not 
include  any  correction  due  to  the  interaction  between  the 
viscous  and  inviscid  regions,  then  the  boundary  layer 
solution  would  be  singular  at  the  point  of  separation. 

In  the  present  work,  three  methods  were  developed  to 
solve  for  subsonic  flows  with  viscous/inviscid  interaction;  a 
finite  difference  method,  an  integral  boundary  layer  method, 
and  a  Hilbert  integral  method.  The  presence  of  the  boundary 
layer  is  assumed  to  affect  the  solution  of  the  inviscid  flow 
region. 
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In  the  first  method,  the  finite  difference  form  of  the 
ANS  equations  are  solved  in  the  viscous  region  while 
Laplace's  equation  written  in  terms  of  the  stream  function  is 
solved  in  the  inviscid  region.  A  coupling  scheme  is 
developed  to  find  the  stream  function  at  the  interface.  One 
cycle  yields  a  new  stream  function  distribution  as  opposed  to 
reference  14 J,  where  many  cycles  are  needed.  Each  cycle  is 
equivalent  to  repeated  Successive  Line  Over-relaxation  (SLOR) 
sweeps.  Edwards  and  Carter  [4]  used  repeated  SLOR  sweeps  for 
each  new  displacement  thickness  to  determine  the  solution  in 
the  inviscid  region.  The  critical  issues  of  this  method  are 
the  efficiency  of  the  viscous  solver  and  the  development  of 
the  fully  implicit  coupling. 

The  second  method  uses  a  finite  difference  scheme  for 
the  inviscid  region  and  an  integral  approach  for  the  boundary 
layer  equations  in  the  viscous  region.  The  displacement 
thickness  now  represents  the  shape  of  a  displaced  body  over 
which  the  flow  is  inviscid.  A  shear  transformation  is 
performed  in  the  inviscid  region  to  allow  for  a  uniform  grid 
in  the  computational  domain.  The  inviscid  solution  will 
produce  and  .  The  Integral  Boundary  Layer  (IBL) 
equations  are  written  in  terms  of  V^,  resulting  in  two 
equations  and  two  unknowns  (5  and  C^)  when  a  velocity  profile 
is  assumed.  The  displacement  thickness  can  now  be  found  and 
used  as  the  next  boundary  condition  to  the  inviscid  solver. 
This  cycle  continues  until  convergence  is  achieved. 
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The  ultimate  goal  of  the  IBL  method  is  the  development 
of  an  integrated  solver  for  the  entire  flow  field.  The  key 
to  finding  an  integrated  solver  is  to  find  an  efficient 
method  for  defining  the  displaced  body. 

The  third  method  solves  the  finite  difference  form  of 

the  boundary  layer  equations  using  the  inverse  mode.  The 

streamwise  velocity  at  the  boundary  layer  edge,  ,  is 

written  as  the  sum  of  the  inviscid  velocity  plus  a  term  which 

accounts  for  the  viscous  effects.  The  viscous  region  is 

solved  to  obtain  U  .  , .  The  Hilbert  integral  finds  the 

e ,  bl 

correction  to  based  on  the  displacement  thickness  used  as 
the  input  to  the  boundary  layer  equations.  The  resulting 
velocity,  will  be  compared  with  If  and 

bl  equal  within  a  set  tolerance,  then  the  method 

continues  until  convergence  is  achieved. 

The  model  problem  chosen  to  develop  the  current  methods 
is  Howarth  115J  flow  over  a  flat  plate,  which  prescribes  a 
piecewise  linear  external  velocity  profile  as  shown  in 
Figure  l.  The  variation  of  the  external  velocity  with  x 
implies  a  pressure  gradient  in  the  streamwise  direction. 

Flow  separation  can  occur  if  the  corner  position  Xq  is  chosen 
correctly.  Briley  (16J  solved  this  flow  using  the  full 
Navier-Stokes  equations,  achieving  separation  with  Xq=0.2  for 
Rg=20800.  The  current  work  will  also  use  Xq=0.2  and  Rg=20800. 
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The  flow  region,  depicted  in  Figure  2,  consists  of  a 
two-dimensional  rectangular  area  where  the  flow  is  subsonic 
everywhere  in  the  field.  The  upstream  boundary  begins  at  a 
non-dimensional  length  of  x=0.05  along  the  flat  plate  where  x 
IS  in  the  streamwise  direction  and  Y  is  in  the  direction 
normal  to  the  plate.  The  downstream  boundary  is  at  x=0.489 
after  Briley.  The  outer  boundary  used  in  the  present  work 
was  Y  =  7.23  where  Briley  used  Y  =  5.4  as  the  outer 
boundary.  The  difference  in  the  two  outer  boundaries  is 
present  to  allow  for  adequate  grid  resolution  in  the  current 
inviscid  solver  while  also  assuring  that  the  interface  is 
positioned  above  the  boundary  layer. 

Results  of  the  three  methods  are  presented  and  compared 
to  existing  work.  Solution  accuracy,  memory  requirements, 
and  computation  time  are  discussed  and  recommendations  for 
further  study  are  given. 
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II  Analysis 


Approximate  Navier-Stokes  (ANS)  Equations 

The  development  of  the  ANS  and  IBL  equations  will  start 
from  the  full  Navier-Stokes  (NS)  equations.  For  2-D 
incompressible,  steady  flow  in  the  absence  of  body  forces 
[1^]  the  NS  equations  in  dimensional  form  simplify  to 

u^  +  Vy  =  0  (2.1) 

*  ™y  =  *  ■■  <“xx  *  V’ 

“''x  *  ''''y  ^  "Py^^  *  “  '''xx  *  ''yy* 

Where 

u  =  velocity  in  x  (streamwise)  direction 
V  =  velocity  in  y  (normal)  direction 
p  =  pressure 

P  =  fluid  density  (constant) 

i'  =  fluid  kinematic  viscosity  (constant) 

The  subscripts  denote  partial  differentiation  with  respect  to 
that  variable.  The  energy  equation  is  not  considered  here 
since  Eqns  (2.1) -(2. 3)  and  the  energy  equation  are  uncoupled  for 
incompressible  flow.  Eqs  (2.1)-(2.3)  can  be  put  in 
nondimensional  form  by  defining  the  quantities 


u' 

=  U/Uj 

v  =  v/y^ 

X' 

=  X/L 

y’  =  y/L 

p' 

(2.4 
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where 

=  free  stream  velocity 
L  =  characteristic  length 

Applying  Eqn  (2.4)  to  Eqns  (2.1)-(2.3)  and  dropping 
the  primes  yields 


U  +  V  = 

X  y 

0 

(2.5) 

uu  +  vu  = 
X  y 

-p  +  R  "^(u  +  u  ) 

*^x  e  '  XX  yy 

(2.6) 

uv  +  vv  = 
X  y 

-p  *  R  ”^(v  +  V  ) 

^y  e  '  XX  yy ' 

(2.7) 

where  the  Reynolds  number  is  defined  by 

(2.8) 

Equations  (2. 5) -(2. 7)  are  the  NS  equations  written  in 
nondimensional  form. 

The  BL  equations  are  obtained  by  considering  the  scaling 
law  of  the  relative  magnitudes  of  velocity  components  in  the 
thin  region.  This  leads  to  an  estimate  of  the  nondimensional 
boundary  layer  thickness 


where  the  symbol  '  implies  order  of  magnitude.  This  is  a 
fundamental  result  of  boundary  layer  analysis.  Notice  that 
for  very  large  R^  the  order  of  magnitude  of  the  boundary 
layer  thickness  is  much  less  than  unity.  Therefore  terms  in 
the  NS  equations  that  are  found  to  be  of  order  of  magnitude  6 
can  be  neglected  compared  to  terms  that  are  of  order  unity. 
First  consider  Eq  (2.5),  the  continuity  equation.  The  u 
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velocity  and  the  streamwise  coordinate  x  are  clearly  of  order 
unity  by  definition.  This  implies  that 

(2.10) 

and  since  the  normal  direction  y  is  of  order  then 


"y  ■  ^ 


2.11 


- 1  2 

Since  is  of  order  ^  from  Eq  (2.9),  the  order  of 


magnitude  of  each  term  in  Eqns  f2.6)  and  (2.1)  can  be 
determined.  The  result  is  that  the  x  momentum  equation  (2.6) 
contains  terms  of  order  unity  except  for  the  term 

.  o 

which  is  of  order  The  y  momentum  equation  (2.7)  is  found 

to  have  all  terms  of  order  6  except  for  the  term 

*^e  '^xx 

3 

which  is  of  order  ^  .  These  two  terms  can  be  neglected 
compared  to  the  relative  magnitude  of  the  remaining  terms. 

The  resulting  set  of  equations  is 

(2.12) 

(2.13) 

(2.14) 


u  +  V  =0 
X  y 


-1 


uu  +  vu  =  -p  +  R  u 

X  y  X  ©  yy 


UVx  +  VVy 


-p  +  R  ^  V 

e  yy 


The  Prandtl  boundary  layer  equations  consist  of  Eqns  (2.12) 
and  (2.13).  Since  all  terms  in  Eq  (2.14)  are  of  order  6,  it 
is  omitted  in  the  classic  boundary  layer  approximation.  The 
ANS  equations  retain  equation  (2.14),  allowing  for  a  change 
in  pressure  normal  to  the  streamwise  direction.  These 
equations  will  now  be  written  in  terms  of  the  stream  function 
t  and  the  vorticity  o. 
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The  stream  function  'f  is  defined  for  2-D,  incompressible  flow 
such  that 

u  =  'J'y  (2.15) 

V  =  (2.16) 

and  the  vorticity,  o,  is  defined  as 

‘^  =  -  Uy  (2.17) 

Notice  that  the  continuity  equation  (2.12)  is  automatically 
satisfied  by  the  definition  of  H*.  The  pressure  can  be 
eliminated  from  Eqns  (2.13)  and  (2.14)  by  taking  the  partial 
derivative  of  Eq  (2.13)  with  respect  to  y  and  subtracting  the 
partial  derivative  of  (2.14)  with  respect  to  x.  Four  of  the 
resulting  terms  are  eliminated  by  the  continuity  equation, 
leaving 

'  (2.18) 


VI  fv  -u  )+v(v  -u  )=R  (v  -u  ) 
XX  xy  xy  yy  e  '  xyy  yyy 

which  can  be  simplified  to 


(2.19) 


by  taking  derivatives  of  Eq  (2.17)  and  substituting  them  into 
Eq  (2.18).  The  velocities  u  and  v  were  also  expressed  in 
terms  of  'f'  from  Eqns  (2.15)  and  (2.16).  The  vorticity  can 
also  be  expressed  in  terms  of  in  Eq  (2.17)  to  obtain 


+  f  +  cj  =  0 

XX 

yy 

Equations  (2.19)  and  (2.20)  are  the  ANS  equations  for  2-D 
incompressible  flow  written  in  terms  of  V  and  <•>.  These 
equations  are  the  model  for  the  viscous  region  in  the  first 
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method,  where  finite  difference  representations  are  used. 
The  inviscid  model  is  obtained  directly  from  Eqn  (2.20), 
where  the  vorticity  is  zero.  The  resulting  equation  is 
simply 


'i'  +  V 

=  0 

XX 

yy 

Integral  Equations 

The  Integral  Boundary  Layer  ( IBL)  equation  can  be 
obtained  by  integrating  the  dimensional  form  of  the  x 
momentum  equation 

\  “yy  (2.22) 

from  y=0  to  y=6.  The  pressure  term  can  be  expressed  in  terms 
of  the  velocity  at  the  edge  of  the  boundary  layer,  U^,  by 
considering  the  boundary  conditions 

“yly=6  -  0 

u  I  r  =  0 
yy'y=5 

Imposing  these  conditions  on  Eq  (2.22)  results  in 

Ue  dUg/dx  =  -p^/p  (2.23) 

Substitute  Eq  (2.23)  into  Eq  (2.22)  and  integrate  over  the 
boundary  layer  to  give 

Jq  ^^^X  '^e  ■  ‘"Jp  ^yy  (2.24) 

Which  simplifies  to 
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J  (uUj^  +  vUy  -  dUg/dx)  dy  =  ( 

where  t  is  the  shear  stress  at  the  wall,  defined  by 
w 

’^w  =  “  “yly^O  ' 

The  normal  velocity  component,  v,  can  be  replaced  in  Eq 
(2.25)  by 


(2.25) 


(2.26) 


=  -J„  “x 


(2.27) 


which  IS  a  result  obtained  by  integrating  the  continuity 
equation.  Substitute  Eq  (2.27)  into  Eq  (2.25)  to  obtain 

J  (uu^  -  dy  -  Ug  dU^/dx)  dy  =  (2.28) 

Integrate  by  parts  to  obtain  the  second  term 

6  y  S  S 

I  '“y  L  “x  “'f  =  “el  “x  ■  J„  ““x 


and  substitute • into  Eq  (2.28)  to  get 

J  (2uu^  -  Vx  '  ""e  =  ’V^ 


(2.30) 


which  after  rearrangement  yields 
6  6 

r  a/ax(u(U  -  u))  dy  +  dU  /dx  j  (U  -  u)  dy  =  r  /p 

Jq  ®  ®  '*0 


(2.31) 


The  displacement  thickness  <5  and  momentum  thickness  a  are 


defined  as 


U3  ^ 


=  J  (Ug  -  u)  dy 


(2.32) 


e  =  r  u(u_  - 


u)  dy 


(2.33) 


Substitute  Eqns  (2.32)  and  (2.33)  into  (2.31)  to  obtain 

d/dx  (U  ^f*)  +  <5*  U  dU  /dx  =  i/P  (2.34) 

C  w  c  w 

Now  cast  Eq  (2.34)  into  nondimensional  form  by  introducing 
the  skin  friction  coefficient  C^,  which  has  the  definition 

C,  =  ^  ^w  (2.35) 

L  ^ 

PU/ 

Expanding  the  first  term  of  Eq  (2.34)  and  writing  the  entire 
equation  in  nondimensional  form  gives  the  final  form  of  the 
integral  boundary  layer  equation 

(2.36) 

An  alternative  form  to  Eqn  (2.36)  is  found  by 

2 

multiplying  it  by  U_  and  rearranging  to  get 

(2.37) 

Eqn  (2.37)  is  used  in  the  second  method.  The  unknowns  are 

C^,  e,  and  5*.  By  assuming  a  velocity  profile,  the  unknowns 

can  be  reduced  to  two  (C^  and  5).  The  second  equation  is  a 

form  of  the  continuity  equation 

u  +  V  =0 
X  y 

Add  and  subtract  the  term  dU^/dx  and  rearrange  to  get 

V  =  dU  /dx  “  dU  /dx  -u  (2.38) 

y  6  c  X 
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Integrating  this  equation  from  y=0  to  y=5  and  using  Eqn 
(2.32)  results  in 

(2.39) 

Equations  (2.37)  and  (2.39)  are  solved  simultaneously. 

Representation  of  the  invlscid  flow  using  a  panel  method 

The  Hilbert  Integral  represents  the  correction  to  the 

inviscid  velocity  due  to  the  effects  of  the  displacement 

thicl^ness.  To  find  an  expression  for  this  correction  Su,  one 

can  consider  a  line  of  source  distribution  q(x)  at  y=0  that 

is  constructed  such  that  the  resulting  flow  ta)ces  place  over 

★ 

a  displaced  body  defined  by  the  displacement  thic)tness  o  (x). 
Figure  3  shows  a  typical  source  distribution  that  results  in 

W 

flow  over  the  body  defined  by  o  (x).  The  source  strength  is 
defined  by  the  constant  mass  flow  rate  it  generates.  By 
considering  a  circular  control  volume  centered  about  the 
source,  the  strength  can  be  shown  to  be  related  to  the  flow 
velocity  at  a  distance  r  from  the  source  by 

q(^)  dF,  =  2nr  V  (2.40) 

The  X  component  of  velocity  d(5u)  for  a  source  of  strength 
q(^)  dl  located  at  x=5  (Fig.  3)  is  simply 

d(5u)  =  V  cos  e  =  q(^)  d5  (X  -  €)  /  2nr^  (2.41) 

where 

r^  =  (x  -  5)^  +  5*^  (2.42) 
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(2.43) 


The  total  5u  considering  the  range  between  and  is 
<5u(x)  =  [  "  q(F. )  (X  -  5)  /  2nr‘'  dE 

Integrating  Eqn  (2.43)  from  zero  to  y  will  give  the  stream 
function  This  integration  results  in 

tan'^(y/(x-E) )  /  2Jr  ]  dE  +  C  (2.44) 

* 

The  relation  between  at  y=5  and  at  y=0  can  be 
approximated  by  considering  a  Taylor  series  expansion  of  M', 
where 

'('(x,0)  ss  '♦'(6*)  -  (**)  =  -  Ug  5*  (2.45) 

since  '^y|y-5*  =  '^(^*)  =  0*  Differentiating  Eqn  (2.45) 

with  respect  to  x  gives 

d'('(x,0)/dx  =  -d(  U.  •''•*) /dx  (2.45a) 

Q 

Differentiating  Eqn  (2.44)  with  respect  to  x  and  taking  the 
limit  as  y  goes  to  zero  and  equating  the  result  to  Eqn 
(2.45a)  will  give  the  simple  relation 

q(.E)  =  2  d(U^^*)/dE  (2.46) 

Substituting  Eqn  (2.46)  into  Eqn  (2.43)  gives  the  final 
Hilbert  integral  expression  for  5u  as 

(2.47) 


The  streamwise  velocity  is  calculated  from 

U  j^(x)  =  PY(x)  +  6u(x)  (2.48) 

which  Includes  the  correction  5u  due  to  the  viscous  effects. 
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Boundary  Conditions 

The  boundary  conditions  will  now  be  discussed  in  general 
for  the  three  methods.  Figure  4  shows  the  computational 
domain  and  governing  equations  for  the  viscous  and  inviscid 
regions.  The  domain  consists  of  a  2-D  rectangular  region 
with  subsonic  flow  throughout.  Boundary  conditions  for  and 
w  are  specified  at  each  of  the  boundaries  (a)-(e). 

At  the  outer  flow  boundary  (a)  the  velocity  is  specified 
by  Howarth's  profile. 

PY(x)  =  1  -  X  ,  for  X  <  Xp 

PY!x)  =  Xq  ,  for  X  i  Xq  (2.49) 

where  X^  is  chosen  to  be  0.2  after  Briley  [16]. 

The  inlet  flow  (b)  is  assumed  known  from  boundary  layer 
theory.  Howarth  [15J  solved  this  flow  using  series 
representations  to  obtain  solutions  up  to  the  separation 
point.  The  resulting  streamwise  velocity  distribution  is 
compared  in  Figure  5  to  the  distribution  obtained  from 
boundary  layer  code  [18J.  The  variable  ETA  in  Figure  5  is 
defined  as  n,  where 

n  =  0.5  y  (2.50) 

Note  the  excellent  agreement  between  the  data.  The  velocity 
obtained  from  the  boundary  layer  code  can  be  integrated  with 
respect  to  the  normal  direction  to  obtain  Y  at  the  inlet. 
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The  inlet  vorticity  is  obtained  from  Eqn  (2.20)  assuming 


=0  in  the  boundary  layer,  resulting  in 


w| 

'x=xl 


-H*  I 

yy ' x=xl 


The  no  slip  condition  prevails  at  (c),  where 

''*y=0  =  \ly=0  =  ° 


and 


vl  =  -'f'  I 
'y=0  x'y=0 


and  therefore 


'^1  =  0 
'  y=0 


Also,  since  't'=o  everywhere  along  the  plate  surface 

I  „  =  0 
xx' y=0 


(2.51) 


(2.52) 


(2.53) 


(2.54 


(2.55) 


(2.56) 


and  therefore  the  vorticity  is  given  by 

'y=0  yy'y=0 
from  equation  (2,20). 

At  the  downstream  location  (d)  a  boundary  condition  on  '1' 
alone  is  required  since  the  ANS  equations  contain  the 


elliptic  term 


xx 


If  the  full  Navier-Sto)tes  equations  were 


used,  the  term  w  would  be  present  in  Eqn  (2.19),  requiring 

AX 

a  downstream  boundary  condition  for  the  vorticity  in  the 

viscous  region.  The  boundary  condition  on  V  is  obtained 

assuming  that  (d)  is  far  downstream  from  the  pressure 

gradient  disturbance.  Uniform  flow  is  assumed  where 

vl  I  n  =  0  (2.57) 

'x=x2  x'x=x2 
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At  the  interface  (e)  the  stream  function  is  obtained 
directly  from  the  coupling  scheme  for  the  finite  difference 
method.  The  integral  and  Hilbert  methods  assume  that  the 
interface  represents  the  displaced  body  over  which  the  flow 
is  inviscid.  For  this  case  'i'=0  along  the  interface.  The 

vorticity  is  assumed  to  be  zero  along  the  interface.  This 
represents  a  departure  from  normal  boundary  layer  type 
conditions  where  the  free  stream  velocity  is  allowed  to  be 
satisfied  asymptotically  at  an  infinite  distance  normal  to 
the  wall  [16].  In  the  present  work,  these  conditions  are 
imposed  at  a  finite  distance  from  the  body  surface. 

The  placement  of  the  interface  from  the  wall  can  be 

quite  important.  Briley's  solution  to  the  present  case  of 

* 

Xp=0.2  shows  that  the  maximum  value  of  o  is  about  0.015. 
For  the  simple  case  of  a  linear  velocity  distribution,  the 
boundary  layer  thickness  is  equal  to  twice  the  displacement 
thickness  [17].  Thus,  for  the  present  work  the  location  of 
the  interface  should  be  at  or  greater  than  approximately 
y=.03,  which  corresponds  to  Y  =  4.32  for  R  =  20800. 


Ill  Methods  of  Solution 


Method  I 


Solution  Procedure.  The  governing  equations  for  this 
method  are  the  ANS  equations  (2.19)  and  (2.20)  for  the 
VISCOUS  region  and  Laplace's  equation  (2.21)  in  the  inviscid 
region.  The  flow  is  initialized  everywhere  by  the  inlet 
conditions.  The  finite  difference  form  of  Eqns  (2.20)  and 
(2.19),  respectively,  can  be  written  as 


A,  +  C,  .H'.  .  ,  +  D,  .iJ.  .  ,  +  E-  tj.  . 

ID  l,D-l  1]  1,D  ID  l/D+l  ID  l,D-l  iD  1,D 


3.1) 


and 


^2d  i,D-1  2d  i,D  2d  i,D+1  2d  i,D-1  2d  i,D 


+  =  RHS^  . 

2d  1,D+1  2,D 


(3.2) 


Where  the  indices  i,D  represent  the  computational  domain  grid 

points  in  the  streamwise  and  normal  directions  respectively. 

The  finite  difference  form  for  Eqn  (2.21)  can  be  written  as 

A-,  .  +  B-.'J'.  +  C-.’f'.  .  ,  =  RHS_  .  (3.3) 

3d  1,D-1  3d  1,D  3d  1,D+1  3,D 

The  calculation  of  the  above  coefficients  and  right  hand 

sides  is  the  subject  of  the  next  two  sections.  The  results 

can  be  found  in  Eqns  ( 3 . 18 ) - ( 3 . 21 )  for  Eqn  (3.3).  The 

coefficients  for  Eqn  (3.1)  will  be  similar  to  the 

coefficients  for  equation  (3.3)  since  the  two  equations 

differ  only  by  the  vorticity  term. 
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The  coefficients  for  Eqn  (3.2)  can  be  found  in  Eqns 

( 3 . 44 ) - ( 3 . 50 ) .  For  now,  only  the  method  of  solution  is 

discussed. 

Successive  Line  Over-Relaxation  (SLOR)  is  used  in  both 
regions.  The  solution  procedure,  depicted  in  Figure  6, 
starts  at  the  inlet  and  marches  downstream.  At  each  i 
location  the  coefficients  of  Eqns  (3.1)-(3.3)  are  found.  For 
the  inviscid  region  this  results  in  a  tridiagonal  system  of 
equations  which  can  be  solved  using  a  form  of  the  Thomas 
algorithm  [19]  described  in  Appendix  A.  This  algorithm  puts 
a  tridiagonal  matrix  of  equations  in  upper  triangular  form 
anc  then  performs  back  substitution  to  find  the  solution 
directly.  From  this  algorithm,  the  solution  at  the  interface 
J  can  be  written  in  terms  of  the  solution  at  the  J+1  point  by 
the  relation 


\1/  •  fi  +  V 

i,J+l  ■  i,J+l  i,J  *i,J+l 


Where  and 


(3.4) 

are  recurrence  coefficients.  The 
boundary  conditions  are  enforced  by  the  definition  of  the 
coefficients  A^ ,  ,  and  RHS^  at  the  boundaries.  For 

the  viscous  region  Eqns  (3.1)  and  (3.2)  result  in  a  block  2x2 
set  of  equations  that  is  solved  using  a  special  form  of  the 
Thomas  algorithm  for  two  partial  differential  equations. 

This  algorithm  is  also  described  in  Appendix  A. 
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The  solution  for  the  viscous  equations  will  have  the  form 


'F 


“i.j  = 


+  s. 

+  T, 

1,3 

1,3  +  1 

1,3 

1,3  +  1 

1,3 

+  s  * 

o 

+  T., 

2,3 

1,3  +  1 

•-  f  J 

1 ,3  +  1 

2,3 

where  R, 


T-,  .  are  the  recurrence  coefficients 

/  J 


'1,  j  ■ 

Eqns  (3. 4) -(3. 6)  the  solution  to  both  the  viscous  and 


(3.5) 

(3.6) 
With 


inviscid  regions  can  be  found  given  all  of  the  boundary 
conditions.  A  crucial  element  of  the  solution  process, 
however,  is  the  coupling  scheme  between  the  two  regions.  The 
governing  equations  at  the  interface  is  given  by  Eqn  (2.21), 
since  the  vorticity  is  assumed  to  be  zero.  The  finite 
difference  equations  at  the  interface  j=J  has  the  form 

AH'  ,  ,  +  bH'  ,  +  cH',  ,  ,  +  dH'. 

i~l,J  i+l,j  i,J+l 


" "  h,j-i 


0  (3.7) 

Where  the  coefficients  A-E  are  known.  The  goal  is  to  be  able 
to  solve  for  H’  in  Eqn  (3.7).  The  value  of  H*  is  known 

1/J 

because  the  solution  is  sweep  downstream.  The  value  of 
H'  is  approximated  to  the  value  from  the  previous 

X  ^  1  f 

level.  Note  that  this  assumption  is  valid  since  the 

difference  in  H'  between  successive  levels  goes  to  zero  as  the 

solution  converges.  The  remaining  unknowns  are  H'  S' 

and  H*  Using  Eqns  (3.4)-(3.6)  the  terms  and 

H'  _  -  can  be  written  in  terms  of  H'  from  the  recurrence. 

1  ^  J  •  1  X  # 

Since  the  vorticity  is  assumed  to  be  zero  at  the  interface, 
the  recurrence  relation  for  given  by  Eqn  (3.5)  can  be 
reduced  at  the  interface  to 
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\U  -  R  V  +  T 

i,J-l  ~  1,J-1  i,J  1,J-1 


(3  .R) 


Substituting  Eqns  (3.4)  and  (3.8)  into  Eqn  (3.7)  results  in 
one  equation  with  one  unknown.  The  stream  function  at  the 
interface,  is  calculated  and  used  as  the  new  boundary 

1  ,  J 

condition.  The  solution  in  the  inviscid  region  is  calculated 
directly  since  it  uses  a  non-recurrence  form  of  the  Thomas 
algorithm.  See  Appendix  A.  Eqns  (3.5)  and  (3.6)  are  used  to 
generate  the  solution  in  the  viscous  region  once  the 
recurrence  coefficients  are  known.  Both  solutions  use  't'  at 
the  interface  calculated  from  Eqn  (3.7)  as  a  boundary 
condition.  The  algorithm  moves  downstream  to  the  next  i 
location  where  the  process  repeats  until  i  reaches  IMAX-1 . 
Residuals  are  then  calculated  and  compared  to  a  tolerance 
value.  If  the  residual  tolerance  is  not  satisfied  the 
iterations  continue  back  at  i=2  and  the  solution  process 
continues  until  convergence  is  achieved. 

Inviscid  Region.  The  model  for  the  inviscid  region  is 
Laplace's  equation 

(3.9) 

Appendix  B  contains  all  of  the  finite  difference  forms  of  the 
derivatives;  allowing  for  a  non-uniform  grid.  The  coordinates 
(x,y)  represent  the  physical  domain  and  (C,n)  represent  the 
uniform  computational  domain.  The  indices  of  the  uniform 
grid  are  denoted  by  i  in  the  5  direction  and  j  in  the  n 
direction. 


'¥  +4*  =0 

XX  yy 
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Eqn  (3.9)  written  in  finite  difference  form  becomes 

. 

D  i,D+l  D  i.3  3  1,3-1  3  1+1,3 

+  E  *  'I'  ,  =0 

3  1-1,3 


where , 


N’  =  "y,:  ^ 


=  -  "y,:  '  "y. 3*1/2  *  \.3-i/2'^ 
^x,j  ^  ^x, 3  +  1/2  ■"  '"'x, 3-1/2^/ 


(3.10) 


(3.11 


(3.12 


"3 

'^y,  3 

3-1/2 

/ 

in'- 

(3.13) 

* 

“3 

=  ^x.3 

^X, 3+1/2 

/ 

45' 

(3.14) 

•* 

^3 

=  'x,3 

^x,:-l/2 

/ 

(3.15) 

as  given  in  Appendix  B.  It  is  implied  that  the  stream 
function  m  Eqn  (3.10)  is  written  at  the  unknown  level 
denoted  as  n+1,  whereas  n  represents  the  current 
level.  The  n+1  level  is  implied  throughout  this  report. 
Terms  at  the  current  level  will  carry  the  n  superscript. 
Rearrange  Eqn  (3.10)  to  obtain 


+  C  *  'F.  .  ,  +  d/  'F. 


'*’i,3  "  *3  '^i,3  +  l  '  ''3  ‘i,!-!  '  *'3  ’i+l,! 


*  =3  ’'1-1.3^  /  ®3 


(3.16) 


Relaxation  is  now  introduced  to  enhance  the  convergence  of 
the  solution.  The  relaxation  parameter  W  is  defined  such 
that  ^  at  the  unknown  level  is  written  as 

(3.17) 


n 


’l,3  =  I  1  -  «  )  ’1,3"  *  "  1,3 
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Notice  that  when  W  =  1,  Eqn  (3.17)  becomes  a  simple  identity 
with  no  relaxation  of  the  solution.  Over-relaxation  (W  >  1) 
can  be  used  to  accelerate  the  convergence  if  the  solution  is 
relatively  stable.  Under-relaxation  (W  <  1)  is  used  mainly 
to  maintain  numerical  stability  in  the  iteration  process. 

In  this  case  more  of  the  solution  at  the  current  value  of 
is  used  to  help  stabilize  the  solution.  Substituting  Eqn 
(3.161  into  Eqn  (3.17)  and  rearranging  gives  Eqn  (3.3) 
rewritten  as 


'V  ,  +  +  C-.'F.  =  RHS_  . 

3j  1,3-1  33  1,3  33  1,3+1  3,3 


where , 


'‘33 

=  W  C  .  /  B 

3  3 

(3.18) 

*33 

=  1 

Hr  * 

(3.19) 

■^3  3 

=  ‘"‘l  /  *3 

(3.20) 

""=3,5  =  '  1  -  «  ‘>3*  's’ 

-  W  e/  r.j  j/  Bj*  (3.211 

All  of  the  terms  in  RHS_  are  known.  The  'V.  term  is  known 

^  1.  f  J 

since  it  is  at  the  current  level.  The  stream  function  at  i-1 
is  known  because  the  solution  is  marched  downstream.  The 
stream  function  at  i+1  is  assumed  to  be  equal  to 
an  approximation.  The  boundary  conditions  are  imposed  by 
redefining  the  coefficients  of  Eqn  (3.3)  at  the  boundaries. 
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The  condition  at  the  outer  boundary  is  given  by  Eqn  (2.49), 
rewritten  as 

’'yly=y2  =  <3-221 

where  PY{x)  is  defined  by  Howarth’s  velocity  profile  as 

PY(x)  =  1  -  X  ,  for  X  <  Xq  (3.23) 

PY(x)  =  Xq  ,  for  X  i  Xq  (3.24) 

Writing  Eqn  (3.22)  in  finite  difference  form  gives 

^32  <  '^i,:2  ■  '^i,d2-1 
which  can  be  written  as 

'^i,j2  =  '*'1,32-1  ^  /  '^yD2> 

Substituting  this  equation  into  Eqn  (3.3)  written  at  j=32-l 
gives 

A  H'  4-  R  W  + 

^3,32-1  1,32-2  ^  “3,32-1  i,32-l 

''3,32-1  '  ^,3^-l  *  ''  '’yd2>'  =  ''”=3,32-1 

which  can  be  cast  back  into  the  form  of  Eqn  (3.3).  The 
resulting  coefficients  are  given  in  terms  of  their  previous 
definitions  as 


^3,32-1  =  ^3,32-1 

®3,32-l  =  ®3,32-l  ^3,32-1 

^«S,d2-1  =  *^«S3,32-1  -  ^3,32-1  /  "y32 


^3,32-1  ° 


(3.28) 

where  the  equal  sign  above  implies  replacement  of  the  terms 
as  in  the  FORTRAN  programming  language.  Also  note  that 
C3  j2-i  zero  after  calculating  the  other  terms. 


31 


The  boundary  condition  for  '1'  at  the  interface  is  given  by  the 
solution  to  Eqn  (3.7)  as  previously  discussed.  For  this  case 
the  coefficients  are  found  to  be 


®3,J+1  *  ®3,J+1 
'^3,J  +  1  ~  "3,J+1 

With  all  of  the  coefficients  known,  the  tridiagonal  solver  is 
used  to  generate  the  recurrence  relation  at  the  interface  and 
to  generate  the  solution  for  The  calculations  follow  the 
solution  procedure  discussed  in  the  previous  section. 


Viscous  Region.  The  viscous  region  is  bounded  by  the 
lower  plate  surface  (Y=0),  the  upstream  and  downstream 
boundaries  at  x=.05  and  x=0.489  respectively,  and  the 
interface  at  Y=4.73.  The  ANS  equations  (2.19)  and  (2.20) 
are  the  mathematical  model  for  this  region.  These  equations 
are  scaled  in  the  vorticity  to  allow  for  the  unknowns  and 
to  be  of  similar  order  of  magnitude.  Define  the  vorticity  as 


O  =  U  /  U) 


(3.30) 


where  w  is  the  vorticity  at  the  plate  at  x=0.05.  This 

ir 

scaling  only  affects  Eqn  (2.19)  since  w  will  become  a  common 
factor  in  Eqn  (2.20).  The  remaining  discussion  will  assume 


the  scaled  form  of  the  equations,  written  as 

f  +¥  +(jo)*  =  o 

XX  yy 

f(i)-'Fw-u)  /R=0 

y  X  X  y  yy'  e 


(3.31) 

(3.32) 
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where  the  vorticity  in  (3.31)  is  the  scaled  vorticity.  The 
convention  used  in  the  previous  section  where  no  superscript 
implies  the  unknown  n+1  level  and  the  superscript  n  implies 
the  current  level  is  retained  here.  Note  that  Eqn  (3.32)  is 
nonlinear  since  the  n+1  level  is  implied  on  each  term  in  the 
equation.  Linearizing  Eqn  (3.32)  results  in 


y  X 


X 


X 


-cj  /R  = 


(3.33) 


yy  e  y  x  x  y 

The  right  hand  side  of  Eqn  (3.33)  as  well  the  n  level  terms 

are  known  from  the  current  level.  Finite  difference 

expressions  for  H'  T  u  and  w  are  obtained  from 

X  y  X  y 

Appendix  B.  The  finite  difference  form  for;  Eqn  (3.31)  is 
exactly  the  result  obtained  from  the  previous  section  with 


the  exception  of  the  vorticity  term.  The  resulting  finite 
difference  equation  is  simply  Eqn  (3.1) 


^lj'^i,D-l  ®lj'’'i,3  ‘^I3'’’i,3  +  l  *  ®lj'^i,j-l  ^Ij^^i,: 


(3.1) 


where  the  coefficients  and  are  given  by  the 

right  hand  sides  of  Eqns  ( 3 . 11 ) - ( 3 . 13 )  respectively  and  the 
remaining  coefficients  are 


(3.34) 

'id  = ' 

(3.35) 

Fjj  =  0 

(3.36) 

RHS,  .  =  -D  .  '1' ,,  .  -  Ea  ^ 

l,j  1  1+1,3  1  1-1,3 

(3.37) 
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The  finite  difference  form  for  Eqn  (3.33)  will  now  be 
addressed.  The  solution  will  be  marched  in  the  streamwise 
direction.  Therefore,  and  o  are  written  as  backward 

X  X 

differences.  The  terms  and  are  written  as  central 
differences.  When  is  negative,  the  convective  term 
will  be  written  as  upwind  difference  to  honor  the  local 
streamwise  direction  in  the  reversed  flow  area.  Now  let 

(3.38) 

(3.39) 


UP  =  1/2  (  ) 

UM  =  1/2  (  ) 


and  notice  that  when  >  o,  UP  equals  and  UM  equals  zero. 
When  '^y  <  0/  UP  equals  zero  and  UM  equals  The  forward 

and  backward  finite  difference  forms  of  w  are 

X 


O  1 

=  P, 

( 

o .  .  - 

x'  1 

XI 

1,3 

o  1 

=  e  . 

{ 

CJ  .  , 

x‘  1 

XI 

1  +  1,3 

• 

and  therefore  the  convective  term  can  be  written  as 

=UP^.  (cj.  -o.,  .  )/A^ 

y  X  XI  '  1,3  1-1,3 


(3.40) 

(3.41) 


.  UM 


OJ  .  ,  .  -  w ,  .  )  / 

1+1,3  1,3 


(3.42) 

Using  Appendix  B,  Eqn  (3.33)  in  finite  difference  form  is 
«xi  '  “l,d  -  "i-1,3  )  /  *  UM  5^1  (  -  “1,3  )  / 

*  “x"  ‘  ”13  "i,3^1 
-  V  '  ”l3  “l.3« 


23  1,3 

-  P  'F  ) 

33  i,3-l  ’ 

P-. 

u 

-  p  u  ) 

33  i,3-l  ^ 

23 

1,3 

'F  ) 

i-1,3  ' 

/  ] 

-  r  Q,  j  -  Q-.  w-  •  +  Qo  •  a  .1  1  /  R 

*■  Ij  i,3  +  l  23  1,3  33  i,3-l  ■'  « 

y  X  X  y 


(3.43) 
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The  terms  ...  are  given  in  Appendix  B.  Equation 

(3.43)  IS  now  cast  into  the  form  of  Eqn  (3.2),  repeated  here 
as 

'f .  .  ,  +  B.,  +  D.,  u  ,  +  E.,  .w. 

23  i,3-l  23  1,3  23  i,3+l  23  1,3-1  23  1,3 


*  f2A.3»i  = 

by  rearrangement  of  the  terms  in  Eqn  (3.43).  The 
coefficients  A-  . . .  RHS_  .  are  given  by 

^  J  t  J 

A  =  -  (1)  ^  p 

23  X  ^33 


(3.2) 


®23 

=  u  ^ 

X 

^23 

y 

^xi  / 

*^23 

=  u  " 

X 

^3 

°23 

X 

**33 

Q33 

/  ^e 

^23 

=  'xl 

(  UP 

-  UM 

)  /  - 

^23 

-  ^ 
“  "  X 

^3 

■  ®13 

/  % 

X  *^23  "  "^23 


U  =  ?  .  (  UP  u ,  ,  .  -  UM  ,  .  )  / 

2,3  XI  '  1-1,3  1+1,3 

y  XI  1-1,3  y  X  X  y 


(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 

(3.49) 

(3.50) 


Eqns  (3.1)  and  (3.2)  are  solved  using  a  special  form  of  the 
Thomas  algorithm  for  two  partial  differential  equations. 

This  algorithm,  which  is  described  in  Appendix  A,  generates 
the  solution  to  a  block  2x2  set  of  equations  in  terms  of 
recurrence  coefficients.  For  the  viscous  region,  the 
solution  will  have  the  form 

...  +  +  T,^  (3.51) 


’'j  ^  *^lj  '’’i,j  +  l  ^13  “i,j  +  l  "^13 
"j  =  ^23  '’'i,j  +  l  ^23  “i,j+l  "^23 


(3.52) 
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The  terms  R.  ...  T-.  are  the  recurrence  coefficients. 

13  j 

Boundary  conditions  are  written  in  the  form  of  Eqns  (3.51) 
and  (3.52)  to  find  the  recurrence.  With  the  recurrence  )cnown 
at  a  boundary,  the  block  tridiagonal  algorithm  can  generate 
the  remaining  recurrence  coefficients.  There  are  two  ways  of 
solving  the  problem.  The  first  way  uses  the  wall  boundary 
conditions  to  generate  the  recurrence  and  the  interface 
conditions  to  generate  the  solution  from  Eqns  (3.51)  and 
(3.52).  The  second  way  is  to  use  the  boundary  conditions  at 
the  interface  to  generate  the  recurrence  and  use  the  wall 
conditions  to  generate  the  solution.  Since  the  current 
coupling  scheme  finds  the  solution  at  the  interface,  the 
first  approach  is  taken.  The  wall  boundary  conditions  are 
given  by  Eqns  (2.54)  and  (2.56).  Since  is  zero  at  the  wall 
and  and  o  at  j=2  are  non-zero,  in  general,  it  follows  from 
Eqn  (3.51)  that 

=  0  (3.53) 

S^^  =  0  (3.54) 

Tj^  =  0  (3.55) 

The  other  recurrence  at  the  wall  are  found  by  expanding  Eqn 
(2.56) 

‘^1  «  =  -'*'  I  n 

y=0  yy'y=0 

into  the  form  of  Eqn  (3.52)  and  comparing  the  coefficients  of 
'V  and  u. 
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The  resulting  recurrence  are 

R22  =  -<  Q31  Qii  >  /  (3.56) 

S22  =  0  (3.57) 

T22  =  0  (3.58) 

With  the  recurrence  known  at  the  wall,  the  block  tridiagonal 
solver  can  generate  the  remaining  recurrence,  where  each 
block  is  a  2x2  matrix.  Eqns  (3.51)  and  (3.52)  are  then  used 
to  find  the  solution  from  the  boundary  condition  at  the 
interface . 

Method  II 

Solution  Procedure.  The  governing  equations  for  this 
method  are  Eqns  (2.37)  and  (2.39)  in  the  viscous  region  and 
Laplace's  equation  (2.21)  in  the  inviscid  region.  The 
inviscid  solver  is  similar  to  the  inviscid  solver  of  method  I 
with  one  exception.  The  lower  boundary  for  this  method 
takes  the  shape  of  the  displacement  thickness  o  as  shown  in 
Figure  7.  The  inviscid  solver  now  solves  for  flow  over  a 
displaced  body  where  '('  =  0 .  The  remaining  boundary 
conditions  are  unchanged.  A  shear  transformation  is  applied 
to  Laplace’s  equation  to  allow  for  a  uniform  grid  in  the 
computational  domain.  The  solution  procedures  is  as  follows: 

1)  Assume  an  initial  5*.  solve  the  inviscid  region 
using  the  SLOR  algorithm. 

2)  Evaluate  and  along  the  displaced  body  from  the 
inviscid  solution. 
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3)  Use  a  Runge-Kutta  method  120J  to  solve  for  5  and 
in  Eqn  (2.37)  and  Eqn  (2.39),  repeated  here  as 

d/dx[  ]  +  Ug  (du^/dx)  5*  =  c^/2 

V  =  d/dx  [  U  6*  ]  -  6  dU^/dx 

* 

The  momentum  thickness  0,  and  the  displacement  thickness  5 

are  written  in  terms  of  5  and  from  Eqns  (2.32)  and  (2.33) 

assuming  a  streamwise  velocity  profile  based  on  the  shear 

stress  at  the  wall.  This  results  in  two  equations  and  two 

unknowns.  Note  that  U  ,  V  and  dU^/dx  are  known  from  Step  2. 

e  e  e 

* 

4)  Repeat  steps  (1)  through  (3)  until  6  does  not 


change . 

Inviscid  Region.  The  model  for  the  inviscid  region  is 

W  +  Y  =0 
XX  yy 

A  shear  transformation  given  by 
5  =  X 

n  =  y  -  6  /  H(x) 


H(x)  =  y^  -  5 


(3.59) 

(3.60) 

(3.61) 

is  applied  to  allow  for  a  uniform  grid  in  the  computational 
domain.  The  transformed  Laplace  equation  becomes 

where 


a=2H^(  1-n)  /H 

A 

^  (  1  -  )  +  i]  / 

C  =  {  1  -  n  )  (  H  -  2 


(3.63) 

(3.64) 

(3.65) 
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Finite  difference  forms  of  the  derivatives  in  Eqn  (3.62)  are 


given  as 


'J'cc  =  (  -  2  'V .  +  '9  ,  .  )  / 

1+1,3  1,3  1-1,3 


(3.66) 


3  'P  ,  -  4  H'  .  ,  ,  +  , 

1,3+1  1-1, 3+1  1-2, 3+1 


/  4  An  A? 


-  [  2  'V  -  4  H*  ,  + 

1,3-1  1-1, 3-1 


+  -i  1  )  /  4  An  A5 

A  •  A.  ,  J  “  -1. 


'f 


nn 

t  =  (  H* 


-  2  t  +  'F  .  ,  )  /  An' 
1,3+1  1,3  1,3-1 


■  ,  )  /  2  An 

1,3+1  1,3-1 


(3.67) 

(3.68) 

(3.69) 


where  a  three  point  backward  difference  scheme  is  used  for 


when  1  >  2,  A  two  point  backward  difference  scheme  is 


used  at  i  =  2.  Substituting  Eqns  ( 3 . 66 ) - ( 3 . 69 )  into  Eqn  (3,62) 


and  rearranging  gives 


N*  ^o-^  *  ®4*  *  =4*  =  '"‘=4’ 


(3.70) 


Where 


-3  u  ^  _§_ 


4  An  A^  An  2  An 


-2 


2  ^ 


B, 


A?‘ 


An‘ 


^  J-  ^ 

,2 


4  An  A^  An  2  An 


(3.71) 

(3.72) 

(3.73) 


'«®4*  =  -  "1-1, 3^"^'  *  “ 


(3.74) 

A  similar  form  of  the  equations  are  obtained  for  the  case  of 


Relaxation  is  now  added  in  the  same  manner  as  the  finite 
difference  inviscid  solver.  Solve  for  'f'.  in  Eqn  (3.70)  and 

1  f  J 

substitute  into  the  right  hand  side  of  the  relaxation 
relation 

'I'  =  (  1  -  W  )  'I'  ”  +  w  *{'  (3.75) 

f  J  /  J  J*  /  J 

to  obtain 


where 


A.  =  W  A.  /  B. 
4  4  4 


^  1 

C4  =  W  c//  b/ 


RHS^  =  (  W  RHS^  /  ®4  )  ^  >  '’'i,j 


n 


(3.76) 

(3.77) 

(3.78) 

(3.79) 

(3.80) 


Boundary  conditions  are  now  applied  to  find  the  values  of 
A^...RHS^  for  the  special  cases  of  j=2  and  j=JMAX-l.  At  j=2, 


N'j  =  2  =  ° 


'y  =  0,  and  therefore 

(3.81) 

At  the  outer  boundary  the  conditions  are  given  by  Howarth's 
velocity  profile,  rewritten  as 


Equation  (3.82)  written  in  finite  difference  form  in  the 
computational  domain  yields 


•'l.JMAX  '  ’’i.JMAX-l  *  «  '''> 

Equation  (3.76)  written  at  j=JMAX-l  is  simply 


(3.83) 


A-  Twikv  o  +  Bil  TWHV  1  +  C-  Tuav  =  *^5.  (3.84) 

4  i,JMAX-2  4  1,JMAX-1  4  1,JMAX  4 
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Substituting  Eqn  (3.83)  into  Eqn  (3.84)  and  rearranging  gives 


the  coefficients  A. ...RHS.  at  j=JMAX-l  as 

4  4 

RHS^  =  RHS^  -  PY(I)  H  An  (3.85) 

B4  =  (3.86) 

=  0  (3.87) 

where  the  order  of  these  calculations  must  be  honored.  With 
all  the  coefficients  known,  the  system  solver  used  in  the 
other  inviscid  solver  can  be  used.  The  solution  procedure  is 
identical  at  this  point  to  the  previous  inviscid  algorithm. 

Viscous  Region.  The  solution  to  the  IBL  equations 
(2.37)  and  (2.39)  start  by  assuming  a  streamwise  velocity 
profile  written  in  terms  of  the  wall  shear  stress.  The 
general  form  is  given  by 

u  =  AU  +  BU  y  +  CU  y^  +  DU  y^  -t-  EU  y^  +  FU  y^  (3.88) 

where  AU-FU  are  determined  by  the  boundary  conditions 


^|y=0  =  ° 

(2) 

(3)  =  0 

( 4  )  U  I  c  =  0 

'  yy|y=5 

151  '‘y|y=0  =  V" 

<5*  “yyy|y=o  =  “ 

A  fifth  degree  polynomial  was  chosen 
a  fourth  degree  form  did  not  provide 


(3.89) 

after  it  was  found  that 
adequate  accuracy. 
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Imposing  Eqn  (3.89)  on  Eqn  (3.88)  results  in 
AU  =  0 
BU  =  r^/iJ 

CU  =  (10/3)  6"^  -  2(t  /n)  6"^ 

DU  =  0 

EU  =  2(t  //J)  6"^  -  5  5"^ 

w  e 

FU  =  (8/3)  -  (r^//J)  5"^  (3.90) 

Eqn  (3.90)  can  be  put  in  non-dimensional  form  using  Eqn 

(2.35),  rewritten  here  as 

2 

C,  =  21-  /Pu 
f  w' 

and  by  considering  the  non-dimensional  forms  given  in  Eqn 
(3.4).  Eqn  (3.90)  in  non-dimensional  form  is  therefore 
AU  =  0 

BU  =  Rg  Cj/2 

CU  =  (10/3)  Ug  5"^  -  5“^ 

DU  =  0 

EU  =  C,  5"^  -  5  U^  5"^ 
e  f  e 

FU  =  (8/3)  Ug  6"^  -  (Rg/2)  5"^  (3.91) 

Eqn  (3.91)  substituted  into  Eqn  (3.88)  defines  the 
non-dimensional  form  of  the  streamwise  velocity  u.  Eqn 
(3.88)  can  now  be  substituted  into  Eqns  (2.32)  and  (2.33)  to 
obtain  expressions  for  the  momentum  and  displacement 
thicknesses  in  terms  of  and  6. 
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While  these  calculations  for  and  6  are  straightforward, 
they  are  also  very  long  and  tedious  for  the  case  of  finding 
the  momentum  thickness  6  (See  Appendix  C) .  The  results  are 


given  as 

6*  =  A  <5  -  B  (3.92) 

where , 

A  =  4/9  (3.93) 

B  =  Rg/30  (3.94) 

and 

0=c6+DC,  6^u'^-E  5^  u  (3.95) 

t  e  re 

where , 

C  =  0.115440115  (3.96) 

D  =  8.297258e-03  (3.97) 

E  =  1.695526e-03  R^^  (3.98) 


Eqns  (3.92)  and  (3.95)  can  now  be  substituted  into  the  IBL 
equations 

d/dx[  ]  +  Ug  (dUg/dx)  5*  =  c^/2  (2.37) 

=  d/dx  [  5*  ]  -  5  dU  /dx  (2.39) 

to  obtain  two  equations  in  two  unknowns,  5  and  C^.  Once  the 

* 

equations  are  solved,  5  can  be  found  from  Eqn  (3.92)  to 
obtain  the  new  shape  of  the  displaced  body  for  the  next  cycle 
of  the  inviscid  solver  discussed  in  the  previous  section. 


44 


Two  methods  for  solving  the  IBL  equations  are  Newton's 


method  for  finding  roots  of  algebraic  equations  and  the 
fourth  order  Runge-Kutta  method  for  integrating  first  order, 
ordinary  differential  equations.  Both  were  implemented  for 
the  IBL  method,  but  only  the  Runge-Kutta  method  120J  was  used 
in  the  final  results.  The  IBL  equations  are  written  as 


d^^/dx  =  F^(Cj,5)  (3.99) 

d(C^)/dx  =  (3.100) 

where  F^  and  F.,  are  functions  of  and  6.  The  integration 
starts  at  x.  =  .05  and  ends  at  x.,  =  .489.  The  expressions 
for  F^  and  F2  will  now  be  found.  Eqn  (2.37)  written  in  terms 
of  and  6  becomes 
■.  cu/a^  . 

"  Vex''* 


-  ^'’ex'^f*  =  “ 


(3.101) 


While  Eqn  (2.39)  becomes 

''e  =  ''Ve  *  ''*“ex  '  2BCf56^  -  -  u^^a  (3.102) 

F^  is  found  by  solving  for  in  Eqn  (3.101)  and 

substituting  into  Eqn  (3.102).  This  results  in  an  equation 
of  the  form  of  Eqn  (3.99),  where 

Fjicj.a)  =  [  -Aau^^  -  -  bdCj6\^z-i 

+  (l/2)BCj6^Ug^z”^  -  ABUgS^Ug^z"^  + 

*“ex*  *  ''e  ]  '  ^ 
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and. 


Z  =  D  5^  -  2  E  Cj  5^  (3.104) 

K  =  AU  -2BC,'5  +  BCS^U  ^z~^  +  2BDC.5^U  z"^ 
e  r  e  re 

-  3BECj^6^z“^  (3.105) 

is  found  by  solving  for  6^  in  Eqn  (3.102)  and 
substituting  into  Eqn  (3.101).  The  result  is 

J  '2U^Ug^C5  -  -  2DUgCj5Q 

+  M  (3.106) 

where , 

P  =  A  -  2  B  5  (3.107) 

Q  =  p-1  (  Vg  -  A  6  +  5  )  (3.108) 

M  =  P”^  (  CU  ^b5^  +  2BDU  0^:6^  -3EBC,^6^  ) 

6  8  t  I 

+DUg6^  -  2ECj5^  (3.109) 

With  Fj  and  F2  determined,  Eqns  (3.99)  and  (3.100)  are  solved 

using  a  fourth  order  Runge-Kutta  algorithm.  The  solution 

produces  and  6  for  each  streamwise  location.  The 

* 

displacement  thickness  5  is  calculated  and  used  as  the  input 

to  the  inviscid  solver,  which  produces  a  new  U  ,  U  and  V 

e  ex  e 

from  the  inviscid  solution.  This  cycle  continues  until  the 

» 

convergence  criteria  on  5  is  satisfied. 
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Method  III 


This  method  solves  the  Boundary  Layer  (BL)  equations  in 
the  viscous  region  and  uses  the  Hilbert  integral,  Eqn  (2.47), 
to  obtain  the  viscous  correction  to  the  inviscid  velocity. 

The  ANS  equations  used  in  the  finite  difference  method  are 
reduced  to  the  BL  equations  by  setting  =  o.  The  solution 
procedure  is  as  follows: 

★ 

(1)  Assume  the  displaced  body  defined  by  5  (x). 

(2)  Solve  the  EL  equations  by  solving  Eqn  (2.19)  and 
(2.20)  with  'V  =  0.  The  boundary  condition  at  the  outer 

ic 

boundary  is  known  from  prescribing  6  in  { 1 ) .  Calculate 
Ue  the  streamwise  velocity  at  the  BL  edge. 


(3)  Solve  the  Hilbert  integral 

5u  =  [d(U^5*)/de]  (X-f,)/r^  dP,  (2.47) 


using  trapezoidal  integration.  On  the  first  pass  in  Eqn 

(2.47)  is  set  to  Howarth's  velocity,  PY.  On  later  passes 

is  taken  as  the  current  value,  which  will  be  denoted  as 

Calculate  the  inviscid  BL  velocity  as 

U  .  =  PY  +  6u  (3.110) 

e,h 


(4)  Generate  a  new  6  (x)  based  on  and  using 


the  following  relation 


n+1  . 


=  6 


n 


<  «e,bl  /  «e,h) 


(3.111) 
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Hr 

This  simple  method  for  updating  5  was  used  by  Carter  and  is 

described  in  Reference  [1].  It  was  noticed  that  small 

deviations  in  the  local  tends  to  preserve  the  volume  flow 

* 

rate  per  unit  width  in  the  BL.  Therefore,  o  a  constant. 

A  local  decrease  in  (adverse  pressure  gradient)  causes  an 

Hr 

increase  in  5  and  vice  versa.  Adding  over-relaxation  to  Eqn 
(3.111)  results  in 

5*  n+1^5*  n[^  ^  (Ue,bl/^e,h^  "  (3.112) 

(5)  With  a  new  5*,  repeat  (2)-(3)  until  h 

within  a  set  tolerance  level. 


IV  Results  and  Discussion 


The  three  methods  discussed  in  Chapter  3  were  applied  to 
Howarth's  flow  over  a  flat  plate.  The  flow  is  incompressible 
with  a  Reynolds  number  of  Rg=20800.  Results  are  given  and 
comparisons  made  where  applicable.  Discussion  on  the 
significant  aspects  of  each  method  are  also  given. 

Method  I 

The  flow  geometry  for  this  method  was  shown  in  Figure 
2  of  Chapter  2.  Briley  [16j  solved  this  problem  with  an 
Alternating-Direction-Implicit  (ADI)  scheme  using  the  full 
Navier-Stokes  equations.  The  fundamental  difference  between 
Briley's  approach  and  the  current  approach  is  the  treatment 
of  the  solution  domain.  Briley  solved  the  Navier-Stokes 
equations  within  a  single  region  where  no  distinction  is  made 
between  inviscid  and  viscous  flows.  The  current  solver 
breaks  up  the  solution  domain  into  a  viscous  region  and  an 
inviscid  region  that  are  implicitly  coupled  together.  The 
advantage  of  the  current  method  is  clear  if  the  inviscid 
region  is  much  larger  than  the  viscous  region.  For  this 
case,  the  relatively  simple  inviscid  model  is  solved  over 
most  of  the  domain  while  the  ANS  equations  are  solved  over  a 
small  region  of  that  domain.  Using  the  full  Navier-Stokes 
equations  in  such  a  domain  would  be  quite  expensive  in  terms 
of  computation  time.  However,  Briley  selected  an  outer 
boundary  that  is  very  close  to  the  edge  of  the  boundary 
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layer.  Therefore,  a  comparison  of  computation  times  between 
these  two  methods  will  not  demonstrate  the  utility  of  the 
current  method.  It  can  be  said,  however,  that  the  larger  the 
inviscid  region,  the  greater  the  amount  of  computational 
savings  using  the  current  method  versus  using  the 
Navier-Stokes  equations  in  a  single  solution  domain. 

Briley's  work  was  used  as  a  checkcase  for  the  current 

fc 

finite  difference  method.  The  displacement  thickness  6  from 
Briley  was  prescribed  at  the  edge  of  the  current  viscous 
solver.  The  resulting  coefficient  of  friction  compared  very 
well  to  Briley's  as  shown  in  Figure  8.  The  grids  used  for 
this  comparison  were  (35x30)  for  Briley  and  (151x74)  for  the 
current  viscous  solver.  The  comparison  shows  that  the 
viscous  solver  is  working  properly.  The  entire  current 
method  was  then  run  for  the  case  of  a  coarse  (30x35)  grid 
used  by  Briley,  an  intermediate  (76x74)  grid,  and  a  fine 
(151x74)  grid.  The  inviscid  grid  was  the  same  as  the  viscous 
grid  in  the  number  of  x  points  with  4  grid  points  used  in  the 
y  direction.  All  the  grids  were  uniform.  The  resulting 
and  ->  are  shown  in  Figures  9  and  10  respectively.  The 
difference  in  the  separation  region  is  due  to  the  treatment 
of  the  boundary  conditions.  With  the  current  zonal  technique 
the  boundary  conditions  are  not  imposed  the  same  way  as  in 
Briley's  case. 
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Figure  9.  Skin  friction  coefficient  for  various  grids 


I 

I  Briley  [ISJ  prescribed  zero  vorticity  and  Howarth's  velocity 

I  at  the  outer  edge  of  his  solution  domain.  In  Method  I,  the 

condition  of  zero  vorticity  is  imposed  at  the  edge  of  the 
I  boundary  layer  and  Howarth's  velocity  condition  is  imposed  at 

the  cuter  boundary  of  the  inviscid  region.  It  is  this 

I 

I  distinction  in  the  boundary  conditions  that  accounts  for  the 

* 

difference  in  the  two  solutions.  When  prescribing  5  at  the 
edge  of  the  viscous  region,  as  was  done  for  the  case  shown  in 
I  Figure  8,  the  separation  region  grows  as  the  grid  is  refined. 

However,  as  Figure  9  shows,  the  current  method  produces  a 
I  smaller  separation  region  as  the  grid  is  refined.  This  event 

IS  explained  by  examining  Figure  10,  which  shows  that 

★ 

refining  the  grid  in  the  current  method  produces  a  lower  5  . 

This  results  in  a  lower  as  well.  In  Briley’s  case, 

* 

however,  the  same  <3  is  prescribed  for  each  refinement  of  the 
grid.  The  two  approaches  are  fundamentally  different. 

Streamline  contours  for  the  current  method  are  shown  in 

I 

'  Figure  11  for  the  intermediate  (76x74)  grid.  The  x  axis  was 

scaled  for  plotting  purposes  only,  and  represents  the  plate 
for  x^=0.05  to  X2=0.489.  The  important  result  from  Figure  11 
is  that  the  streamlines  at  the  interface  (7=4.73)  are 
completely  continuous.  The  coupling  scheme  produces  a  smooth, 
continuous  solution  between  the  viscous  and  inviscid  regions. 
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Figure  12  compares  the  predicted  inviscid  velocity  at  the 
interface  using  Method  I  to  Howarth's  velocity.  The 
difference  in  these  two  curves  is  due  to  viscous  effects. 

The  current  method  was  also  run  for  a  higher  value  of 
the  corner  velocity  Xq ,  with  the  expectation  that  a  higher  Xq 
will  result  in  a  longer  region  of  adverse  pressure,  which 
should  result  in  a  larger  separation  region.  This  was 
exactly  the  case  as  shown  in  Figure  13.  Here  X^  was 
increased  from  0.2  to  0.21.  The  resulting  increase  in  the 
separation  region  is  apparent  from  the  figure. 

The  computation  times  for  the  current  method  for  the 
coarse,  intermediate,  and  fine  grids  are  55,  579,  and  4101 
CPU  seconds  respectively  on  the  ASD  CYBER.  The  convergence 
history  for  the  intermediate  grid  is  shown  in  Figure  14.  The 
residual  tolerance  was  lo”^  for  all  grids.  The  approximate 
memory  use  is  14600  words  for  the  intermediate  grid. 
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Figure  lA.  Convergence  history  for  Method  1 
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Method  II 


The  IBL  method  was  described  in  Chapter  3.  The  fifth 
degree  velocity  profile,  along  with  the  IBL  equations,  modeled 
the  viscous  region.  The  inviscid  region  solved  the  finite 
difference  form  of  Laplace's  equation  for  a  non-uniform  grid. 
The  inviscid  grid  is  (76x51),  with  the  outer  boundary  at  the 
same  location  as  the  finite  difference  method.  The  residual 
tolerance  value  is  lO”^.  Relaxation  was  used  in  the  inviscid 
region  where  W  was  finally  chosen  to  be  1.3. 

ic 

The  IBL  method  was  first  executed  with  5  and  U  from 

e 

the  finite  difference  method  used  as  the  initial  conditions. 
For  this  special  case,  the  inviscid  solver  is  not  required 

since  U_  is  given.  The  normal  velocity  can  be  calculated 

♦ 

Knowing  <5  and  U^.  The  resulting  for  only  one  iteration 
of  the  method  is  shown  in  Figure  15.  The  excellent  agreement 
between  this  and  Briley's  demonstrates  that  the  method 
can  be  a  useful  alternative  to  solving  the  viscous  region. 

The  CPU  time  required  to  solve  the  viscous  region  is  0.122 
seconds  on  the  ASD  CYBER.  In  addition.  Figure  15  shows  that 
this  rather  simple  method  is  clearly  capable  of  solving  for 
flows  with  an  adverse  pressure  gradient;  including  flow 
separation.  Further  iterations  of  the  current  IBL  algorithm 
results  in  an  instability  in  the  solution,  which  result  in 
very  inaccurate  values  of  C^  past  x  0,3.  The  primary 
reason  for  this  behavior  is  thought  to  be  the  assumption  of 
the  velocity  profile. 
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Figure  15.  Skin  friction  for  Method  II 
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Figure  16  shows  the  resulting  for  the  case  of  a  fourth 

degree  and  a  fifth  degree  velocity  profile  with  the  same 
» 

inputs  of  5  and  .  Notice  that  the  improvement  in  the 

solution  between  the  two  cases  is  due  to  the  accuracy  of  the 

assumed  velocity  profile.  The  fourth  degree  profile  was 

generated  with  the  same  boundary  conditions  as  the  fifth 

degree  profile,,  except  that  ^yyjy=^  ^  ®  fourth  degree 

profile.  Further  accuracy  in  the  method  should  be  attainable 

by  specifying  a  sixth  degree  profile,  where  '•'yyy|y-6  ~  ® 

would  provide  the  additional  boundary  condition. 

An  additional  run  with  the  IBL  method  was  made  to 

investigate  the  small  ridge  in  the  of  Figure  14  at  x=0.2. 

This  ridge  is  believed  to  be  caused  by  a  slight  irregularity 
* 

in  the  input  <5  .  To  confirm  this  hypothesis,  a  sixth  degree 

polynomial  was  generated  to  have  the  approximate  shape  of  the 

displacement  thickness.  This  6  is  shown  in  Figure  17.  The 

was  obtained  using  the  inviscid  solver.  The  resulting  C^, 

shown  in  Figure  18,  shows  that  smooth,  continuous  output  is 

* 

obtained  from  the  method  when  a  smooth  5  is  used  as  the 

input.  This  verifies  that  the  small  ridge  in  the  of 

Figure  15  is  due  to  an  irregularity  in  the  displacement 

thickness.  This  is  an  important  point,  since  the  inviscid 

* 

solver  uses  first  and  second  derivatives  of  o  to  obtain  the 
transformation  parameters  of  the  computational  grid. 
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Figure  17.  Polynomial  displacement  thickness 
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Figure  18.  Skin  friction  for  polynomial  displacement 
thickness 
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The  total  CPU  time  for  one  cycle  of  the  IBL  method  is 
38.622  seconds,  where  38.5  of  this  time  is  due  to  the 
inviscid  solver  and  the  calling  program.  This  time  can  be 
greatly  reduced  by  using  less  grid  points,  by  reducing  the 
residual  tolerance  to  lO”^,  and  possibly  by  using  an  ADI 
method.  The  memory  requirements  are  approximately  17000 
words,  of  which  11,628  words  are  used  to  store  the 
transformation  parameters  ot,  P,  and  C  over  the  (76x51)  grid. 
These  parameters  are  stored  to  calculate  the  residuals.  A 
scheme  could  be  developed  to  calculate  the  residuals  locally 
in  the  flow  by  saving  only  three  local  columns  of  the  data. 
This  would  reduce  the  total  memory  requirements  down  to 
approximately  5500  words. 

Method  III 

Method  III  was  explained  in  Chapter  3.  The  boundary 
layer  equations  are  solved  over  the  viscous  region.  The 
Hilbert  Integral,  Eqn  (2.47),  is  then  solved  to  obtain  the 
correction  to  the  Inviscid  velocity  due  to  viscous  effects. 
The  method  iterates  until  from  the  boundary  layer 
calculation  matches  obtained  from  the  Hilbert  integral 
calculation.  The  resulting  C^  for  this  method  is  shown  in 
Figure  19  for  two  values  of  Xq.  The  data  near  the  inlet  at 
x^  and  at  the  corner  position  required  smoothing  due  to  a 

weak  instability  at  these  points. 
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Figure  19.  Skin  frict 


This  instability  was  also  noticed  by  Cebeci  and  Stewartson 

121],  who  added  several  modifications  to  their  algorithm  to 

smooth  the  introduction  of  the  Hilbert  integral  near  the 

inlet.  The  most  severe  modification  was  the  addition  of  an 

artificial  correction  term  to  the  Hilbert  integral  that  added 

a  maximum  5u  at  the  inlet  and  added  less  5u  as  x  increased. 

Cebeci  also  multiplied  by  the  term  cosec (( x-x^ ) /2Ax)  in 

the  integration  over  the  first  four  points  to  further  smooth 

the  data.  The  only  modification  to  the  current  method  was  to 

smooth  the  data  near  the  inlet  and  near  the  corner  position 

Xq .  Adding  an  artificial  term  to  6u  may  help  smooth  the 

resulting  data,  however  it  also  directly  adds  an  error  to  the 

solution.  Indeed,  Cebeci 's  result  for  C^  showed  no 

separation  at  all  for  the  case  of  Xq  =  0.21.  Figure  19 

clearly  shows  separation  for  this  case  and  is  in  excellent 

qualitative  agreement  with  Briley's  result.  The  slight 

irregularity  of  the  solution  in  Figure  19  near  the  inlet  is 

of  little  concern.  Additional  smoothing  can  be  done  but 

it  will  not  affect  the  overall  trends  in  the  solution.  The 

displacement  thicknesses  for  the  current  method  are  shown  in 

* 

Figure  20.  They  are  in  qualitative  agreement  with  the  o 
from  the  finite  difference  method  and  from  Briley’s  results. 
The  resulting  boundary  layer  edge  velocities  are  shown  in 
Figure  21.  UEBL  is  from  the  boundary  layer  analysis  and 
UEH  is  from  the  Hilbert  integral  calculation.  As  expected, 
the  two  profiles  match  for  the  converged  solution. 
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Figure  20.  Displacement  thickness  for  Method  III 
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The  CPU  time  for  this  method  was  60  seconds.  An 

over-relaxation  of  the  solution  was  used  where  W=1.3.  The 

convergence  criteria  specified  that  the  maximum  difference 

-4 

between  UEBL  and  UEH  at  any  x  location  be  less  than  15x10 
This  criteria  resulted  in  39  iterations  to  achieve 
convergence.  The  convergence  history  is  shown  in  Figure  22. 
The  ordinate  represents  the  error  in  normalized  from  the 
error  of  the  first  iteration.  The  method  required 
approximately  14900  words  of  memory. 
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Figure  22.  Convergence  history  for  Method  III 
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Comparison  of  Results 


The  CPU  times  and  memory  requirements  are  summarized  In 
Table  1,  Also  presented  are  the  average  absolute  errors 
(taken  over  five  points  in  x)  between  C^  of  the  current 
methods  and  Briley's  C^ .  For  Method  I,  c^  from  the 
intermediate  grid  is  used  in  the  comparison  since  the  other 
methods  have  the  same  number  of  x  points  (See  Figure  9).  The 
skin  friction  coefficient  for  methods  II  and  III  can  be  found 
in  Figures  15  and  19  respectively.  In  Figure  19,  the 


Xq  =  0.2  case  is  used. 

CPU  Time 

Memory 

c 

Mean  f  error 

Method  I 

-course  (30x35) 

55  sec 

-inter.  (76x74) 

579  sec 

14.6  KW 

0.096 

-fine  (151x74) 

4101  sec 

Method  II  (1  cycle) 

-viscous 

0,122  sec 

0 . 1 

-inviscid  (76x51) 

38.500  sec 

5.5*  KW 

Method  III 

60  sec 

14.9  KW 

HP  HP 

0.14 

Table  1.  Summary  data  for  Methods  I, II,  and  III 
*  using  local  scheme  to  calculate  residuals 
**  after  smoothing 
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V  Conclusions  and  Recommendations 


Method  I  uses  finite  difference  equations  with  SLOR 
sweeps  for  solving  the  ANS  equations  in  the  viscous  region 
and  the  stream  function  equation  in  the  inviscid  region.  An 
implicit  coupling  scheme  matches  the  two  solutions.  The 
solutions  obtained  from  Method  I  compared  well  with  full 
Navier-Stokes  solutions.  The  coupling  scheme  developed  for 
this  method  provided  an  efficient  means  of  patching  the 
viscous  and  inviscid  regions.  In  addition,  an  initial 
displacement  thickness  is  not  required  to  start  the  solution. 
The  cycle  time  for  the  course  grid  was  the  fastest  of  any  of 
the  three  methods,  with  the  resulting  skin  friction 
coefficient  very  close  to  that  of  the  finer  grids. 

Method  II  uses  finite  difference  approximations  for 
solving  the  stream  function  equation  in  the  inviscid  region 
and  a  fourth  order  Runge-Kutta  method  for  solving  the 
integral  boundary  layer  equations  in  the  viscous  region. 
Method  II  was  shown  to  give  very  good  for  the  first 
Iteration.  Solutions  for  past  the  first  iteration  become 
less  accurate.  It  is  recommended  that  a  velocity  profile 
specified  by  a  higher  order  polynomial  be  used  to  determine 
if  an  improvement  in  the  stability  of  the  solution  can  be 
achieved.  The  major  contribution  of  this  method  is  the 
efficient  viscous  solver.  The  solution  to  the  boundary  layer 
equations  was  reduced  to  finding  the  solution  to  a  coupled 
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I 

I  set  of  first  order,  ordinary  differential  equations  that  were 

solved  using  a  simple  fourth  order  Runge-Kutta  method. 

Method  III  obtained  the  inviscid  flow  solution  by  a 
panel  method,  while  the  viscous  flow  solution  is  obtained 
using  the  finite  difference  form  of  the  boundary  layer 
equations  operating  in  an  inverse  scheme.  The  solutions 
obtained  from  Method  III  were  in  general  agreement  with  the 
known  solutions. 

For  the  current  model  problem  and  geometry,  Method  I 
provided  the  best  overall  performance  as  evidenced  by  the 
data  given  in  Table  1  of  Chapter  4.  However,  for  more 
complex  geometries,  Method  II  would  have  the  best  potential 
for  providing  efficient  solutions  with  a  minimal  amount  of 
required  memory. 

Any  further  study  of  these  methods  should  consider  the 
extension  to  3-D  flow,  compressible  flow,  and  also  to  flow 
over  more  realistic  geometries. 
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Appendix  A:  Tridiagonal  System  Solvers 


Single  System  of  Equations 

This  algorithm  was  obtained  directly  from  Appendix  A  of 
Reference  [Ij.  It  is  used  in  the  current  work  to  solve  the 
system  of  equations  resulting  from  the  finite  difference  form 
of  Laplace's  equation,  written  in  general  form  as 


A  ,  +  E  't'  c  4'  =  RHS. 

1  l,:-l  D  l-D  D  1,3+1  ; 


al 


where  i  is  fixed  and  3  varies  from  3=2  near  the  interface  to 
3=JM1  near  the  outer  boundary.  Writing  Eqn  (al)  for  all 
values  of  3  results  in  the  matrix  equation 


^,2 

^1,2  •••’ 

0 

»  • 

1,2 

*  • 

RHS., 

\,3 

^.3  •••• 

‘^i,3 

.1,3 

RHS3 

(a2 

0 

^1, JMl 

‘^l,  JM1> 

^  1,  JMlJ 

where  the  terms  A,  «  and  C,  are  written  in  terms  of  the 

1  f  J!  1/UNX 

remaining  coefficients  from  the  boundary  conditions.  The 
algorithm  simply  takes  matrix  (aZ)  and  cast  it  into  upper 
triangular  form  by  performing  a  series  of  row  operations, 
given  as 


®i,3  ^  ®i,3  ■  ^1,3  ^i,3-l^  ®i,3-l 
RHS.  =  RHS.  -  A,^^  RHS,^  ..^/ 


(a3) 

(a4) 


Where  the  equal  sign  implies  replacement  as  in  the  FORTRAN 
programming  language.  The  upper  diagonal  terms  remain 
unchanged  because  of  the  tridiagonal  form  of  the  equations. 
The  equations  can  now  be  solved  by  simple  back  substitution 
starting  at  j=JMl  and  marching  upwards. 
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Block  Tridiagonal  Solver 


This  algorithm  was  obtained  from  Reference  [22J.  The 
general  features  of  the  solver  are  exactly  like  that  of  the 
Thomas  algorithm  described  in  [IJ  except  tor  the  current  case 
a  block  2x2  set  of  equations  exist. 

The  finite  difference  representations  of  the  governing 
equations  can  be  written  in  the  following  general  form 


■  J-1 

w 

D  o 

1  J-1 

F  Cl) 

J-1 

o 

II 

(a5) 

’J-1 

+  B., 

^2 

D  o 

^2  J-1 

+  + 

F  Cl) 

2  J-1 

-  *^2 

(a6) 

where  the  coefficients  are  functions  of  the  grid  and  the  J 
subscripts  denote  the  normal  direction.  Define  the 
recurrence  relations  as 


= 

R 

+  S,  , 

+  T, 

(a7) 

j 

IJ 

J+1 

IJ 

J+1 

IJ 

'^J  = 

^2J 

'V 

J+1 

®2J 

o 

J+1 

'^2J 

(a8) 

Equations 

(a3) 

and 

(a4) 

were 

chosen  this  way  to 

allow  for  a 

marching  scheme  that  calculated  the  recurrence  coefficients 
R.  ,  ,T.,  ,R., ,  S- ,T.,  starting  at  J=0  and  marching  towards 

J=JMAX.  Then  the  solution  'f.,  o  are  found  by  knowing  the 

\J  s) 

boundary  condition  at  JMAX  and  marching  the  solution  towards 
J=0  knowing  the  recurrence  coefficients.  The  marching 
direction  is  arbitrary,  but  choosing  the  inarching  in  this  way 
is  compatible  with  the  development  of  the  viscous/inviscid 
coupling  scheme  used  in  the  present  work. 
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Substitute  ( a7 )  and  (a8)  into  (a5)  and  (a6)  to  get 


-1 

''j 

+ 

=1, 

w  +  T  ) 

J-1  J  1,J-1' 

"  ®1  '^J 

+ 

4' 

J+1 

*'>l'''2,J-lV  ®2,J-1  “j  *  ’'2,J-l' 

*  “j  '  ''l  "j+1  =  =1 

(a9) 

and 

-1 

''j 

+ 

=1, 

J-1  “j  *  '^l,J-l' 

^  "'j 

+ 

tm 

’'j 

.,+D.,(R-  _  .4'  +  S_  _  ,  w  +  T-  -  -  ) 

4 

•  E  u  +  F  W  =  G 

2  J  ^2  J+1  ^2 

(alO) 

Now  rearrange  { 

a9) 

and  (alO)  into  the  form 

+ 

^=1 

4*  +bo+Fw  =d 

J+1  1  J  1  J+1  1 

(all) 

4. 

+ 

^2 

4^  +bo+Fw  =d 

J+1  °2  J  2  J+1  °2 

(312) 

where 

^1 

r 

®1 

4* 

^1 

^^1,J-1  °1  '^2,J-1 

(al3) 

= 

®2 

+ 

^2 

, J“1  ^  ^2  ^2  J“1 

(al4) 

‘^l 

= 

^1 

+ 

^1 

®1,J-1  *  °1  ^2, J-1 

(al5) 

b. 

= 

^2 

+ 

^2 

®1,J-1  °2  ®2,J-1 

(al6) 

^1 

= 

- 

^1 

"^l.J-l  ■  ^1  "^2, J-1 

(al7) 

^^2 

= 

^2 

- 

^2 

'^l.J-l  ■  ^2  "^2, J-1 

(al8) 

Now  multiply  Eqn  (all)  by  b2  and  multiply  Eqn  (al2)  by  b^  and 
subtract.  The  resulting  equation  after  some  rearranging  is 

’’j  =  '"^2  ‘’l  -  ^  'fj.l  *  '<^2  *>1  -  '’l  ‘=2'/“o'  “jtl 

♦  ((dj  bj  •  '*2 

where 

Dq  =  b2  “  ^2  (a20) 
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Compare  Eqn  (al9)  to  Eqn  (a7)  to  find 


^IJ  =  ^^2 

^1 

-  ^1 

(a21 

-  ^ 

(a22 

^2 

-  d_ 

i. 

{a23 

The  remaining  recurrence  coefficients  are  found  by 
multiplying  Eqn  (all)  by  a^  and  multiplying  Eqn  {al2)  by  a^ 
and  subtracting.  The  vorticity  is  rearranged  in  the  form  of 
Eqn  f a8 )  to  find  the  recurrence 


2J  ~  ^^1 

^2  ■  ^2 

a,  : 
1 

(a24) 

2J  = 

a 2  ~  F., 

(a25) 

2J  = 

*^2  ■  ^2 

(a26) 

With  the  recurrence  known,  the  solution  can  be  obtained  from 
Eqns  (a7)  and  (a8).  Notice  from  Eqns  (al3)-  (al8)  that  the 
recurrence  at  the  J-1  level  are  needed  to  generate  the 
solution.  The  specification  of  the  boundary  conditions  at 
the  wall  (J=l)  gives  the  recurrence  for  J=l.  Knowing  this, 
the  remaining  recurrence  are  calculated. 
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Appendix  B:  Finite  Difference  Expressions 


The  general  finite  difference  expressions  for  the 
derivatives  in  the  governing  equations  are  presented  here. 
Because  of  the  inarching  scheme,  the  derivatives  in  the  normal 
fy)  direction  can  be  taken  as  central  difference  while  the 
streamwise  derivatives  are  generally  taken  as  either  forward 
or  backward  difference  depending  on  the  flow  direction.  The 
elliptic  H"  term  is  taken  as  central  difference. 

A  non-uniform  grid  is  assumed  where  the  coordinates 
(x,y)  represent  the  physical  domain  and  (5,n)  represent  the 
computational  domain  with  indices  i,j  respectively.  The 
metrics  of  the  transformation  are  calculated  as 

"yh  =  ‘  "yh-i/s  *  '’yl3-i/2  ’ 


Where , 

'^yh-i/2 

and 


^x|i  ^  ^  ^x|i  +  l/2  ^x|i-l/2^ 

where , 


(b2) 

(b3) 


(b4) 


^x| 1+1/2 
^x|i-l/2 


AF,/(Xi^l  -  X.) 

-  x^_^) 


(b5) 

(b6) 
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The  derivatives  calculated  in  the  physical  domain  can  now  be 
expressed  in  terms  of  the  computational  domain  by  simple 
application  of  the  chain  rule. 


\U  -  <  "yh^i/2  ^  ''yl3-i/2  ' 

=  <  '’ylj-n/2  '*n\i*i/2  *  ^13-1/2  '’'n|3-i/2  • 

=  '’yl3^i/2  '  ^*1  -  "3 

*  M3-1/2  '  ^ 

=  O.W^n  ^3,1  -  l'>y|3.i/2  -  M3-1/2>  "3 

■  '’y|3-3/2  ’'3-1] 


where , 

^3  =  '’■'  M3*:/2  ' 

Pjj  =  -0.5  (  -  Oyij.jyj  )  /  ao 

”33  =  M  3-1/2  ' 

Similarly, 

U  I  =  P,  ,  0)  ,  +  .  U  -  P_  .  (i)  ,  , 

yh  11  1  +  1  2:  j  3d  d-1 

The  first  derivative  expressions  in  x  are  simply 

x|i  ^x|i  '1  1-1  ^ 

u  I  .  —  ^  I  .  ((1),  “(1).  }  / 

x|l  x|l  '  1  i-l 


(b8) 

(b9) 

(blO) 

(bll) 

(bl2) 

(bl3) 
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The  second  derivatives  are  found  as  follows 


T 


yyh  ^  ^yh 


=  n„u  ( 


-  'f 


yh  y|j+i/2  ■  yh-i/2 


)  / 


'^yh^\h+i/2  '^r,\j+i/2  ■  ^yh-l/2  '’'nlD-l/2 

=  Q,  ,  H'  .  -  'f  +  Q,  ,  , 

D+1  ^2j  3  33  3-1 

where , 

«i:  =  "yh  Yh-i/:  ^ 

«2D  =  'Vh  ‘  yh*i/2  *  Mj-in  '  ' 

'h:  =  M:  '■’yh-i/a  ' 

Similarly , 

w  I  =  0,  w .  ,  -  .  o  +  0,  .  w  , 

yyh  Id  d+i  23  3  33  d-i 

The  second  derivatives  in  x  will  have  the  same  form  as 
derivatives  with  the  metrics  interchanged,  therefore 

\x|i  =  ■'li  -  Xai  \  *  '^31  ^-1 


where , 


^li  ”  ^x|i  ^x|i+l/2  • 

*^2i  ^  ^x|i  ^  ’x|i+l/2  ^xli-1/2  ^  ^ 

^3i  ^  ^x|i  ^x|i-l/2  ^ 


)/^n 

(bl4) 

(bl5) 
(bl6  ) 
(bl7) 

(bl8) 
the  y 

(bl9) 

{b20) 

(b21) 

{b22) 
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Appendix  C:  Details  of  Method  II 

it 

The  calculation  of  5  and  ^  in  Eqns  (3.92)  and  (3.95) 
will  be  shown  here  in  greater  detail.  The  development  of  the 
functions  F^(C^,S)  and  F.,(C^,^)  are  also  presented.  Method 
II  assumes  a  velocity  profile  given  by  Eqn  (3.88)  with  the 
non-dimensional  coefficients  given  in  Eqn  (3.91). 

.Substituting  (3.88)  into  Eqn  (2.32),  repeated  here  as 

,<s 

=  I  (l-u/U  !  dy  (cl) 

0  ® 

yields 

'■*  =  L  1  -  0-5F.gCjUg'^y  -(  (10/3)6‘2-R^C^6‘^U^‘^y^ 

-(ReCfUe"^5‘^-55‘^)y^  -  {(8/3)^"^  -  0 . ) y^] 

(C2) 

Integrating  Eqn  (c2)  gives 

a*  =  y  -  0.25R^C£Ug‘^y^-(l/3)[(10/3)5’2-R^C^5‘^Ug*^]y^ 

-(1/5) [R^C^Ug‘^5'3-55'^]y5  ^ 

-(1/6)  [(8/3)6'^-0.5RgCgU^’V^]y®  (C3) 

which  results  in 

6*  =  a  -0.25RgCjUg"^a^-(10/9)5  +(l/3)RgCjUg~^6^ 

-(l/5)RgC^Ug"^62  +  6  -(4/9)5  +  ( 1/12 ) R^C^U^'^a^  (c4) 

Eqn  (c4)  simplifies  to 

6*  =  (4/9)5  -{l/30)Rg 

which  is  Eqn  (3.92)  with  the  terms  A  and  B  from  Eqns  (3.93) 
and  (3.94)  respectively. 
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To  obtain  the  expression  for  B  in  Eqn  (3.95),  substitute  the 
velocity  profile,  Eqn  (3.88),  into  Eqn  (2.33),  repeated  here 
as 


e 


(u/U^)  -(u/Ug)^Jdy 


where  the  velocity  profile  can  be  written  as 
u 

when  BU,  CU,  EU,  and  FU  are  given  by  Eqn  (3.91) 


BU  y  +  CU  y^  +  EU  y^  +  FU  y^ 


{c6) 

(c7  ) 
Squaring 


Eqn  (c7)  results  in 

=  BU^y^  +  2 (BU)CUy^+2 {BU)EUy^  +2 ( ( BU) FU+ ( CU ) EU) y® 

24  7  ">8  9210 

+  CU^y  +  2(CU)FUy  +  EU  y  +  2(EU)FUy  +FU  y 


(c8) 


Substituting  Eqns  (c7)  and  (c8)  into  Eqn  (c6)  results  in 


e  =  J  [  (BU)Ug'^y  +  (cu)Ug‘^y^  +  (EU)Ug‘^y^  +  (FU)Ug‘^y^ 

-(BU)^U^"^y^  -2BU(CU)U^"^y^  -  2BU(EU)U^"^y^ 

Q  ^  cf 

-2(BU(FU)  +  (CU)EU)U^"^y^  -  (CU)^U  ”^y^  -  2CU(FU)U  "^y"^ 

6  6  6 

-(EU)"Ug  y  -  2FU(EU)Ug  y  -  (FU)  U^  -y  J  dy  (c9) 


where , 


(BU)^  =  Re^Cg^/4 

(CIO) 

BU(CU)  =  (5/3)RgCjUg5"^  - 

(l/2)Rg^Cj^6-^ 

(Cll) 

BU(EU)  =  (l/2)Rg^Cj^5"^  - 

(5/2)RgC^Ug6"^ 

(cl2) 

BU(FU)  =  (4/3)R^CjUg6”5  _ 

(l/4)Rg2Cj25-^ 

(cl3) 

CU(EU)  =  (25/3)RgCjUg5“^ 

-(50/3)Ug^6“® 

(C14) 
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BU(FU)+CU(EU)  =  ( 29/3 )RgC^U^5"^  - { 5/4 ) 

-(50/3)U^^5‘®  (cl5) 

(CU)^  =  (  100/9)U^"6"^  -(  20/3)R^CjUg'5‘^ 

{cl6) 

CU(FU)  =  ( 80/9  -  (  13/3)R^CjU^6'‘^ 

+1 1/2>R^^C^^6‘^  (cl7) 

(EU)^  =  R  -lOR  +  25U  ^5"®  (cl8) 

e  t  ere  e 

EU(FU)  =  (31/6)RgC^U^6"®  -  (40/3)Ug^6"® 

-(l/2)Rg^Cj^6’‘^  (Cl9) 

(FU)^  =  (64/9)U^^5‘^°  -  (8/3)RgCjUg5‘® 

+(l/4)Rg^Cj^5"®  (C20) 
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Now  integrate  Eqn  (c9)  and  substitute  in  Eqns  (cl0)-(c20) 


to  obtain 


-1 


B  =  (1/4)  +  (1/3)U^ 

+  (1/5)U  "^[r  C,<5"^  -  5U  5"^]6^ 

A  ^  Cl  T  A  * 


-  R  C.6‘^]' 


+  (  1/6  3/3  )U^cS"^-(  l/2)R^Cj5’^j6'= 

-  (l/3)U^’"[Rg^C^^'4]6^ 

-  ;i/2)U^“^[(5/3',RgC^U^^''^  -  (l/:)R^^Cj-f>‘^]' 

-  (1/3)’J  "^[(l/ZlR^-C.^cS'^  -  (5/2)R^C,U  6"^]' 


e  f 


'e^f^e 

2.-4 


2/7)U^  ‘L(29/3)RgC^U^5  - ( 5 /4 ) '-(50/3)Ug^6  °J6 


2^-4 


-3.„  2„  2, 


-  (1/5)U^  ^1  (  100/9)U^^6  ’-(20/3)R  C^US  "’+R^^C,^6  ^6 


ere 


'e 


-  {  1/4  )U^"^[(80/9)U^^^^"^-(13/3)R^CjUg6‘^+(  1/2)R^^C^^6"^]i 


'e  "'f 


2.-6 


-  ll/9)U^  -lOR^C^U^S  '  +  25U^^6  °Jf- 


-  (l/5)Ug‘"[(  31/6  )RgCjUg5"®-( 40/3  )Ug^5“°-(  1/2  )Rg^C£^6‘'^]6^ 


2.-10 


2„  2, 


-(1/11)U^  ^L(64/9)U^^'?  -(8/3)R^C^Ug^  +  (  1 /4  ) R„  “^'3  '"J5 


e  ''f 


Collect  common  terms  in  Eqn  (c21)  to  obtain 


9  =  L(10/9)-1+(4/9)+(100/21)-{100/45)-{80/36)-(25/9) 
+{8/3)-{64/99)]5  +  [( 1 /4 )-{ 1/3 )+( 1/5 )-( 1/12 )-( 5/6 )+( 5/6 
-  (58/21) +  (4/3 )  +  (13/12)  +  { 10/9) -(31/30)  + (8/33)] 

+[( -1/12 )+( 1/4 )-( 1/6 )+( 10/28 )-( 1/5 )-(l/8)-(l/9)+( 1/10) 

-(1/44)1r.^C^^5®u 


(C23) 


Eqn  (c22)  reduces  to 
0  S!  0.1154401155  +  8.297258e-03 

e  r  e 

-1 .695526e-03  R 

e  r  e 

which  is  simply  Eqn  (3.95)  with  the  values  of  C,D, 
given  by  Eqns  ( 3 . 96 ) - ( 3 . 98 )  respectively. 

The  expressions  for  F^(C^,5)  and  F2(Cj,5)  will 
derived.  The  process  starts  be  writing  Eqns  (2.37) 
(2.39)  in  terms  of  and  5.  This  results  in  Eqns 
and  (3.102),  repeated  here  as 

2U  U  C5  +  CU  ^5  +  DU_.  C.5^  +  DU.S^C,..  +  2DU„C.56 


e  ex 


'ex"f 


fx 


'e"'f  X 


-  ®“ex'^f«  =  “ 


and 

V  =  A'S  U  +  A5u  -  2BCx:55  -  -  U  5 

e  X  e  ex  f  x  fx  ex 

To  find  F^(C^,5),  solve  for  in  Eqn  (3.101)  and 

into  Eqn  (3.102).  The  first  step  results  in 

■^fx  =  t  -''"e'^x  - 

■►3ECj^«^6jj+(l/2)CjUg^-U^U^jjA«+BU^jjCj62]/z 

where , 

2  =  DU  5^  -2EC,5^ 

e  r 


and  E 

now  be 
and 

(3.101) 


(3.101) 

(3.102) 
substitute 

{c24) 

(c25) 
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Substituting  Eqn  (c24)  into  (3.102)  results  in 


V  = 

A6  +  A5U 

-2BC.56 

e 

X  e  ex 

f  X 

-B5^| 

-CU  ^5  z~^  - 

2C6u  U  Z~^  -DC,5^U  z' 

■  e  X 

e  ex  f  ex 

+  3EC^ 

=^^^5xZ‘^+(1/2 

)  CjUg^z"^-UgU^j^A62"^+BU 

f  X  e 


ex  f 


{c26) 


The  term  can  now  be  solved  in  terms  of  and  5  from  Eqn 


( c . 26 ) ,  yielding 


2  2-1 

(C27) 


-BU^AS^z‘^Ug^+B^C^5'*z‘^U^^+U^^6+vJ  /K  (  c27 

where , 

K  =  AUg  -2BC^6+BC<r.^Ug^z"^  +  2BDCj6^UgZ"^-3BEC^^6^z”^  (c28) 

Eqn  (c27)  is  Eqn  (3.103)  of  Chapter  3,  where  5  equals 
Fl(Cj,5)  from  Eqn  (3.99).  is  found  by  solving  for 

in  Eqn  (3.102)  and  substituting  into  Eqn  (3.101).  The 
first  step  results  in 


=  [v^  -  a6u^^  +  6u  +  b5^C.  J/P 

*•6  ex  ex  f  x^ ' 


where , 


P  =  A  U  -  2B  C,6 
e  f 


(c29) 


(c30) 


Eqn  (c29)  can  be  simplified  further  by  grouping  the  terms 


that  do  not  contain  that  is 


"x  =  5  * 


where , 


Q  =  L  -  a6u  +  6u  J  /P 
*•  e  ex  ex^ ' 


(c31) 


(C32) 
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Now  substitute  Eqn  {c31)  into  Eqn  (3.101)  to  obtain 
+  2DUgCja[Q  +  -  2EC^5^  -  3EC^^5^[q  + 

-(l/2)CjUg2+U^Ug^A:S-BU^^C^62  =  0 

Now  solve  for  in  Eqn  (c33)  to  obtain 

^fx  =  [  '^^^^e^ex  ■  ^^e^'2  '  -2DUgC^5Q  +  3ECj^^^ 

+  (1/2)C^U^^  -  BS^C^UgJ/M 

where , 

M  =  ECU  ^6^p”^+DU  '5"  +  2BDU  C,5^p'’^-2EC.6^-3BEC,^5^p"^ 
e  e  e  r  t  r 

Eqn  (c34)  is  Eqn  (3.106)  of  Chapter  3,  where  equals 

F.,{C^,6)  from  Eqn  (3.100). 


] 

(c33) 

■Q 

(c34  ) 

(c35) 
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ABSTRACT 


The  study  of  flows  with  viscous/inviscid  interaction  has 
attracted  many  researchers  over  the  last  decade.  These  flows 
occur  whenever  the  adverse  pressure  gradient  is  large  enough  to 
cause  flow  separation.  The  current  emphasis  is  to  find  efficient 
ways  of  solving  these  types  of  flows  without  solving  the  full 
Navier-Stokes  equations. 

Three  methods  for  solving  the  viscous/inviscid  problem  were 
studied.  The  first  method  uses  finite  difference  equations  to  model 
both  the  viscous  and  inviscid  regions.  A  coupling  scheme  is  developed 
to  match  the  two  solutions.  The  second  method  solves  the  integral 
boundary  layer  equations  in  the  viscous  region  and  finite  difference 
equations  in  the  inviscid  region.  The  third  method  solves  the  Hilbert 
integral  to  generate  a  correction  to  the  inviscid  velocity  using  the 
boundary  layer  equations  as  the  viscous  model.  The  model  problem  used 
in  this  work  is  Howarth  flow  over  a  flat  plate. 

The  three  methods  were  evaluated  in  terms  of  solution  accuracy, 
memory  requirements,  and  computation  times. 


